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Abstract 

Biofluidic  microsystems  are  essentially  a  miniaturized  microchip  implementation  of  an  ana¬ 
lytical  chemistry  laboratory.  Such  microsystems  are  fast,  accurate,  automatable  and  inexpensive  to 
fabricate  and  have  spanned  sample  collection,  preparation  and  analysis  for  biological/chemical 
agent  detection  and  on-chip  drug  or  chemical  synthesis  and  delivery  for  warfighter  defense. 

The  design  of  such  biofluidic  microsystems  has  posed  a  significant  challenge.  An  initial 
generation  of  design  tools  based  on  continuum  simulation  and  reduced  order  modelling  was  devel¬ 
oped  through  DARPA  support  in  the  late  ‘90s.  These  tools  were  useful  for  component  design  and 
optimization,  however,  they  were  too  cumbersome  to  be  applied  at  the  system  level.  As  designers 
evolved  from  designing  individual  channels,  pumps,  valves  and  mixers,  to  whole  assays  on  a  chip, 
this  shortcoming  became  critical. 

This  report  describes  approaches  to  system-level  design  tools  and  methodologies.  It  incorpo¬ 
rates  significant  advances  in  a  behavioral  modeling  and  simulation  methodology  for  Lab-on-a-Chip 
(LoC)  design.  This  modeling  and  simulation  methodology  was  developed  with  optimization  for 
design  synthesis  in  mind,  enabling  rapid  automated  design  of  highly  complex  biofluidic  microsys¬ 
tems. 

The  behavioral  modeling  methodology  involves  decomposing  a  complex  LoC  system  into  a 
small  set  of  elements.  Each  of  the  elements  is  associated  with  a  parameterized  behavioral  model 
that  describes  its  electric  and  biofluidic  behavior.  Key  issues  addressed  include  schematic  represen¬ 
tation,  behavioral  multi-physics  modeling  and  numerical  and  experimental  validation.  The  model¬ 
ing  effort  focused  on  sample  transport  in  LoC  devices.  Turn  and  Joule  heating  induced  dispersion  in 
electrophoretic  separation  chips  are  modeled  using  the  method  of  moments.  The  skew  of  the  species 
band  is  effectively  represented  by  a  set  of  Fourier  cosine  series  coefficients  that  are  obtained  analyt¬ 
ically.  These  coefficients  capture  the  effect  of  band  skew  on  separation  performance  in  various 
complex  chip  geometries  (including  multiple  turns).  Variations  of  sample  concentration  profiles  in 
laminar  diffusion-based  micromixers  are  also  derived  using  the  Fourier  cosine  series  representation. 
The  model  holds  for  arbitrary  sample  flow  ratios  and  inlet  concentration  profiles,  and  accurately 
considers  the  overall  effects  of  device  topology,  size  and  electric  field  on  mixing  performance.  In 
addition,  a  simplified  reaction  model  is  developed  and  integrated  with  the  separation  and  mixing 
models  to  perfonn  system-level  schematic  simulations  of  an  integrated  LoC  system.  Simulation 
results  at  both  element  and  system  levels  are  validated  against  numerical  and  experimental  data. 


ii 


Excellent  accuracy  (generally  less  than  <  5%  in  relative  error)  and  tremendous  speedup  (>  100X) 
have  been  achieved  when  compared  with  finite  element  analysis.  The  resulting  modeling  and  simu¬ 
lation  framework  is  a  significant  contribution  to  balancing  the  needs  for  efficiency  and  accuracy 
thus  enabling  iterative  design  of  complex  biofluidic  LoC  systems. 

These  models  were  integrated  in  a  simulation  backplane  that  captured  the  complex  physio- 
chemical  phenomena  and  enabled  the  development  of  design  tools  that  simultaneously  considered 
these  phenomena  as  well  as  challenging  chip  layout  and  channel  interconnectivity  issues.  The 
chemistry  that  takes  place  during  chip  operation,  as  well  as  the  chip  layout  and  manufacturing  pro¬ 
cess  must  be  understood  so  that  the  appropriate  design  trade-offs  and  constraints  are  considered. 
Channel  geometry,  and  the  system’s  channel  topology  are  shown  to  contribute  a  great  deal  to  the 
overall  perfonnance  of  the  final  LoC  design.  These  issues  result  in  a  design  problem  that  is  highly 
nonlinear  and  highly  combinatorial. 
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(a)  BioMEMS  microneedles  for  transdermal  drug  delivery  [2],  (b)  Multiplex  Lab-on-a-Chip  DNA 
sequencer  integrates  384  independently  operatable  capillary  electrophoretic  separation  channels  for 
high  throughput  analysis  [3],  (c)  Gene  expression  microarray  with  >  1  million  probes  with  zoom  in  of 
individual  probe  locations  after  hybridization  reaction  [4].  5 

(a)  Cross  section  of  EWOD  device  [12];  (b)  Droplet  injection  from  reservoir  (left)  and  split  (right);  (c) 
Independent  and  simultaneous  droplet  transport;  and  (d)  two  droplets  are  merged  (left  and  center),  and 
then  merged  with  reservoir  (right)  [14].  7 

(a)  Spiral  Capillary  Electrophoresis  separation  assay  for  amino  acids  [20];  (b)  Microreactors  for  drug 
discovery  [21];  (c)  Multifunction  immunoassay  that  integrates  mixing,  reaction,  and  separation  [22]. 


System-level  Immunoassay  schematic  showing  straight  channels,  U  turns,  elbow  turns,  and  other  ele¬ 
ments  used  to  integrate  mixing,  reaction,  injection  and  separation  functions  on  a  single  LoC  [41]. 

10 

Related  modeling  methods.  13 

(a)  The  activation,  zk,  of  a  node  is  a  function  of  the  weighted  inputs  of  the  previous  layer's  activations. 

(b)  General  feedforward  neural  network  topology  utilizing  multiple  hidden  layers  with  an  unspecified 

number  of  nodes  with  non-linear  activation  functions,  and  a  single  output  node  using  a  linear  activation 
function.  14 

k-fold  cross  validation  algorithm.  First  divide  the  N  points  in  the  dataset  into  k  subsets  (3  to  5  subsets  is 
typical).  Then,  identify  one  subset  ( i )  as  the  validation  set,  used  to  compare  with  the  trained  neural  net¬ 
work.  The  other  k-1  subsets  (  ~i ),  are  used  to  train  the  neural  network.  Each  subset  is  the  validation  set 
exactly  once.  The  error  for  each  iteration  is  averaged  at  the  end  to  provide  an  estimate  of  the  generaliza¬ 
tion  error  for  the  network.  This  is  an  efficient  use  of  data,  since  all  data  is  used  for  training  and  valida¬ 
tion.  There  is  no  exclusive  validation  set.  16 

(a)  Uniform  sampling  on  a  grid  in  two-dimensions  with  a  total  of  16  data  points  resulting  in  effectively 
4  data  points  per  dimension,  (b)  Uniform  sampling  off-grid  in  two  dimensions,  again  with  16  data 
points,  resulting  in  effectively  16  points  per  dimension,  a  better  more  denser  sampling  per  dimension. 

17 

A  neural  network  is  trained  on  the  function  shown  on  the  left.  The  network  has  a  single  hidden  layer  of 
tansig  neurons  and  one  linear  output  layer.  It  is  trained  for  450  max  iterations,  and  the  sample  size  is 
varied.  The  error  measured  is  the  mean  squared  error  of  the  network  to  the  training  data  (not  a  measure 
of  the  generalization  error).  The  Sobol  quasirandom  sequence  is  the  best  performing,  followed  by  the 
Gaussian  sequence,  and  finally  the  uniform  grid.  The  Gaussian  sequence  is  a  simpler  algorithm  and  is 
compatible  with  domain  shaping  methods.  19 

(a)  Behavioral  model  structure  for  the  converging  intersection  in  the  micromixer.  Index  1,  r  and  out  rep¬ 
resent  the  quantities  at  the  left  inlet,  right  inlet  and  outlet.  Both  electric  (V)  and  biofluidic  (dm)  pins  are 
defined  at  the  terminals  of  the  model.  Electrically,  it  is  modeled  as  a  combination  of  three  resistors  (Rl, 
Rr  and  Rout)  with  zero  resistance.  Different  sample  concentration  profiles,  cl(h)  and  cr(h),  at  inlets  are 
merged  and  compressed  at  the  outlet  cout(h).  (b)  Behavioral  model  structure  for  the  diverging  intersec¬ 
tion  in  the  micromixer.  Similarly,  index  1,  r  and  in  represent  the  quantities  at  the  left  outlet,  right  outlet 
and  inlet.  Sample  concentration  profile  at  inlet,  cin(h)  is  split  and  stretched  out  into  two  parts,  cl(h)  and 
cr(h),  flowing  out  of  the  outlets.  22 

Behavioral  model  structure  for  separation  channels  in  electrophoretic  separation  microchips.  Index  in 
and  out  represents  the  quantities  at  the  inlet  and  outlet.  Both  electric  (V)  and  biofluidic  (t,  Sn  ,  s2  and  A) 
pins  are  defined  at  the  terminals  of  the  model.  Electrically,  the  channel  is  modeled  as  a  resistor  (R).  The 
variations  of  biofluidic  pin  values  due  to  dispersion  effects  are  captured  by  the  model.  25 
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Verilog-A  description  for  a  180x  turn  involving  clockwise  flow  of  the  species  band.  It  determines  the 
signs  used  in  Eq.  (24),  as  well  as  the  canceling  and  strengthening  effects  of  the  skew.  28 

(a)  An  electrokinetic  focusing  micromixer.  Sample  a,  flowing  from  the  top  input  channel  to  the  intersec¬ 
tion,  is  pinched  by  sample  b  (or  buffer)  from  both  side  channels.  Then  samples  mix  in  the  bottom  mix¬ 
ing  channel,  (b)  Its  hierarchical  schematic  representation.  The  triple  input  and  one  output  cross 
intersection  is  modeled  as  a  cascade  connection  of  two  converging  intersections.  30 

(a)  Schematic  simulation  results  (lines)  compared  with  numerical  data  (symbols)  on  concentration  pro¬ 
files  c  (sample  a)  along  the  normalized  channel  width  h  for  the  electrokinetic  focusing  and  T-type  mix¬ 
ers.  In  contrast  to  the  T-type  mixer,  the  focusing  mixer  considerably  improves  sample  homogeneity  due 
to  the  reduced  diffusion  distance  between  samples;  (b)  Schematic  simulation  results  on  variation  of 
mixing  residual  Q  along  axial  channel  length  (data  points  are  connected  by  lines  to  guide  the  eye)  for 
the  electrokinetic  focusing  mixer  involving  different  stream  width  s.  A  smaller  stream  width  (e.g., ) 
yields  a  lower  initial  mixing  residual  and  a  more  uniform  concentration  profile  at  the  end  of  the  mixing 
channel.  30 

An  electrokinetic  serial  mixing  network  [74]  and  its  hierarchical  schematic  representation.  The  network 
consists  of  reservoirs,  mixing  channels,  T-  and  cross-intersections.  Sample  and  buffer  is  released  and 
collected  by  the  reservoirs.  In  the  composable  approach,  the  serial  mixing  network  is  represented  as  a 
collection  of  interconnected  mixing  elements  composed  of  microchannels,  converging  intersections  and 
diverging  intersections.  32 

Comparison  of  experimental  data  [79]  with  DC  schematic  simulation  on  variance  s2  vs.  separation  time 
t  of  species  a  in  a  serpentine  electrophoretic  separation  microchip  of  two  complementary  turns.  The 
grey  bars  represent  the  residence  time  of  the  sample  within  the  turns.  The  first  turn  skews  the  species 
band  and  accordingly  incurs  abrupt  increase  in  variance.  The  transverse  diffusion  in  the  inter-turn 
straight  channel  smears  out  most  of  the  skew  and  presents  a  nearly  uniform  band  before  the  second  turn 
(see  the  inset  of  numerical  simulation  plot).  The  second  turn  then  distorts  the  band  again  in  the  opposite 
direction,  leading  to  another  turn-induced  variance.  34 

Transient  analysis  simulates  the  electropherograms  output  from  three  detectors,  which  are  respectively 
arranged  before  the  first  turn  (top  trace),  in  the  middle  of  the  inter-turn  straight  channel  (middle  trace) 
and  after  the  second  turn  (bottom  trace).  Attributed  to  the  difference  in  electrokinetic  mobility,  the  spac¬ 
ing  between  concentration  peaks  of  species  a  and  b  increases  as  they  migrate  through  channels.  The  dis¬ 
persion  effect  leads  to  the  continuous  decreases  in  the  band  amplitude.  34 

(a)  A  spiral  electrophoretic  separation  microchip  [83].  It  consists  of  five  turns  with  continuously 
decreased  radius  (1.9,  1.8,  1.7,  1.6  and  0.8  cm).  Within  them,  species  Dichlorofluoroscein  flows  in  the 
same  direction  (clockwise),  (b)  Its  hierarchical  schematic  representation.  35 

Comparison  of  schematic  simulation  results  to  the  experimental  data  on  theoretical  plate  number  Ns  vs. 
electric  field  E.  Right  axis  shows  the  relative  error  between  simulation  and  experiments.  The  linear 
growth  of  the  plate  number  Ns  with  electric  field  E  implies  that  molecular  diffusion  is  the  major  disper¬ 
sion  source  in  such  a  system.  36 

(a)  A  hybrid  electrophoretic  separation  microchip.  It  consists  of  both  spiral  and  serpentine  channels. 
Species  flows  in  the  clockwise  direction  in  both  turns  T1  and  T2  (spiral  topology),  thereby  T2  strength¬ 
ens  the  sharp  skew  generated  by  Tl.  The  skew  almost  persists  through  the  inter-turn  straight  channel 
between  T2  and  T3  and  is  significantly  cancelled  out  by  T3  (serpentine  topology),  (b)  Its  hierarchical 
schematic  representation.  36 

Comparison  of  numerical  data  with  DC  schematic  simulation  on  variance  vs.  separation  time  in  the 
hybrid  electrophoretic  separation  microchip.  Very  sharp  skew  (see  Fig.  20)  is  generated  and  the  vari¬ 
ance  accumulates  after  turns  T 1  and  T2  due  to  their  spiral  topology.  The  skew  almost  persists  through 
the  inter-turn  straight  channel  between  T2  and  T3  and  is  significantly  cancelled  out  by  T3  attributed  to 
their  serpentine  topology,  which  as  a  result  yields  a  drastic  variance  drop  after  T3.  37 

Basic  operating  stages  for  various  injector  topologies.  The  most  common  topologies  are  the  cross,  dou¬ 
ble-tee,  and  gated  cross.  The  earliest  on-chip  injectors  were  the  tee.  38 
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Two  species  are  combined  and  a  band  of  material  representative  of  an  injector  formed  plug  is  placed  at 
the  beginning  of  the  channel.  Species  one  and  two  have  a  diffusivity  of  1x10-10  m2/s.  Species  one  has  a 
mobility  of  5.25x10-7  m2/Vs,  whereas  species  two  has  a  mobility  of  5.5x10-7  m2/Vs.  After  traveling  in 
a  50000  V/m  electric  field  for  1.2mm,  the  separation  of  the  two  species  bands  is  measured.  In  the  first 
example  (a),  the  plug  has  a  width  of  lOum.  In  the  second  example  (b),  the  plug  has  a  width  of  50um. 
Even  though  the  plug  in  (b)  contains  more  material,  it  is  not  as  easy  to  resolve  as  the  smaller  plug  in  (a). 

39 

Proper  operation  of  an  injector  requires  the  identification  of  control  parameters  and  their  applicable 
ranges.  Here  identifying  the  fields  (arrows)  is  the  first  step,  but  then  the  magnitude  and  ratios  must  be 
additionally  constrained  for  proper  operation  (  E1L  >  0,  E1D/E2D  >  a  ,  where  E1D=E3D)  41 

Leakage  results  for  the  cross  and  double-tee  injectors.  The  values  of  the  loading  stage  parameters  are 
fixed,  but  different  in  each  of  the  four  plots.  The  blue  region  defines  infeasible  operation,  and  the  red 
region  defines  feasible  operation  for  the  dispensing  stage  field  ratios,  eD  and  Peclet  number  PeD. 

43 

Comparison  of  actual  injection  plug  to  its  effective  Gaussian  representation.  The  actual  output  of  a  dou¬ 
ble-tee  injector  is  in  the  top  channels,  the  effective  Gaussian  model  is  shown  in  the  bottom  channels. 
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1.  Summary 

One  of  the  original  DARPA  supported  biofluidic  microsystem  researchers,  Jed  Harrison,  of 
the  University  of  Alberta  said  that  “It  will  become  important  to  place  design  capability  within  the 
hands  of  the  chemist  or  biologist  who  understands  the  application  if  we  are  to  see  a  dynamic  leap 
forward  in  the  acceptance  of  micro  fluidic  technology”  [1].  This  has  guided  project  proposal  as  well 
as  its  implementation. 

The  primary  goals  of  the  “Synthesis  of  Biofluidic  Microsystems”  project  were  to  develop: 

•  Shorten  the  development  cycles  of  biofluidic  microsystems;  at  the  time  of  the  proposal  quoted 
times  for  custom  biofluidic  chip  design  was  on  the  order  of  1  year.  The  proposal  goal  was  to  cut 
it  down  to  a  few  days,  or  a  greater  than  10X  improvement. 

•  Enable  design  of  more  complex  biofluidic  microsystems;  at  the  time  of  the  proposal,  the  state  of 
the  art  was  a  single  separation  system  or  a  single  reaction  system  with  a  few  channels  with  one 
or  two  biomolecular  species.  The  proposal  goal  was  to  enable  multiple  parallel  assays  with 
more  than  10  biomolecular  species  and  to  integrate  thousands  of  channels,  or  a  greater  than 
100X  improvement. 

•  Develop  a  lab-on-a-chip  design  methodology  that  was  primarily  focused  on  synthesis  to  impact 
the  development  of  structured  iterative  design  methodologies. 

All  the  goals  were  met  and  exceeded  during  the  course  of  the  project.  For  example,  the  tools 

and  methodologies  developed  in  this  project  were  used  to  design  a  chip  with  more  than  4  meters  of 
microchannels  to  simultaneously  perfonn  multiple  assays  that  ranged  from  separating  DNA  mole¬ 
cules  and  separating  amino-acids  to  an  immunoassay  that  combined  both  separation  and  reaction 
subsystems  in  less  than  one  hour! 

To  achieve  these  innovations  we  have  developed  a  modeling  and  simulation  methodology 
that  partitions  complex  lab-on-a-chip  systems  into  commonly  used  elements.  Parameterized  models 
are  developed  for  each  of  these  elements  that  no  longer  need  artificial  parameter  fitting  from  contin¬ 
uum  simulation  that  used  to  be  the  state  of  the  art  at  the  time  this  project  started.  Also,  we  have 
developed  mathematical  fonnulations  for  optimal  design  of  LoC  subsystems  and  multifunction 
LoCs.  A  rigorous  general  disjunctive  programming  fonnulation  for  the  optimal  design  of  multi¬ 
plexed  LoCs  as  well  as  a  more  tractable  combinatorial  fonnulation  was  also  developed. 

These  LoC  system  design  tools  have  a  direct  insertion  path  in  the  recent  drive  to  low-cost 
sequencing  engines  needed  for  further  advancement  in  genomics  and  proteomics.  The  DoD  shared 
this  interest  for  use  in  DNA  dogtags  for  warfighter  identification.  Additional  applications  are 
expected  for  the  design  of  a  variety  of  point  of  care  devices  for  hand-held  whole  blood  monitors  and 
bio- weapons  detectors. 
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2.  Introduction 

Microdevices  are  finding  increasing  application  in  biology  and  chemistry,  and  the  life  sci¬ 
ences.  A  driving  biomedical  instrumentation  vision  involves  extracting  a  fluid  sample  from  the 
human  body,  analyzing  it  to  diagnose  disease,  and  then  injecting  a  drug  to  treat  it,  all  within  the 
same  microdevice.  Constituent  components  for  this  vision  have  been  introduced  over  the  last 
decade,  and  can  be  partitioned  into  three  categories.  BioMEMS  devices  are  made  using  planar  IC 
manufacturing  techniques,  augmented  with  MEMS-based  sacrificial  etching  processes.  BioMEMS 
devices  are  typically  discrete,  such  as  valves,  pumps  or  needles  (Fig.  1(a)).  Lab-on-a-Chip  (LoC)  or 
Micro  Total  Analysis  System  (pTAS)  devices  are  biofluidic  integrated  circuits,  with  fluid  instead  of 
current  flows  (e.g.  DNA  separation  assay  in  Fig.  1(b)).  A  third  class  of  device,  called  BioChips  or 
Microarrays  (Fig.  1(c)),  integrates  surface  chemistry  needed  for  biomolecule  identification  (bio¬ 
sensing).  The  discrete  BioMEMS  devices  are  not  really  chips,  and  will  not  be  considered  in  this 
paper.  The  planar  integrated  BioChip  and  Lab-on-a-Chip  platforms  are  both  suitable  for  biosensing, 
and  can  use  the  BioMEMS  devices  to  interface  to  the  human  body  for  sampling  fluid  from  the  body 
and  delivering  a  drug  to  the  body. 

Biosensing,  or  analysis  of  unknown  biological  markers  in  a  sample,  is  a  basic  lab  function. 
Let  us  consider  a  benchtop  lab  to  see  what  needs  to  be  integrated  on  a  chip.  Analysis  begins  with  the 
sample  being  mixed  with  known  fluorescently  labelled  probes  (using  shakers,  stirrers  or  centri¬ 
fuges).  This  mixture  is  then  converted  into  a  reaction  product,  typically  with  the  help  of  catalysts, 
heat,  light,  or  enzyme.  Next,  the  reaction  products  are  placed  in  a  field  (e.g.  electric  or  gravitation) 
separating  the  constituents  according  to  their  charge  or  mass.  Whether  the  reaction  took  place  or 
not  can  thus  be  identified  using  the  fluorescent  marker,  and  leads  to  the  lab  technician  labelling  the 


FIGURE  1.  (a)  BioMEMS  microneedles  for  transdermal  drug  delivery  [2],  (b)  Multiplex  Lab-on-a-Chip 
DNA  sequencer  integrates  384  independently  operatable  capillary  electrophoretic  separation  channels  for 
high  throughput  analysis  |3],  (c)  Gene  expression  microarray  with  >  1  million  probes  with  zoom  in  of 
individual  probe  locations  after  hybridization  reaction  [4]. 
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test  as  positive  or  negative. 

BioChips  use  chemical  reactions  to  identify  unknown  biomolecules.  A  library  of  probes  are 
immobilized  at  specific  locations  on  the  chip  surface  to  construct  the  BioChips.  Fluorescently- 
labelled  unknown  bioagents  react  with  these  probes.  The  amount  of  fluorescence  at  each  location  is 
analyzed  by  image  processing  and  reveals  the  relative  quantitation  of  several  thousand  different 
nucleic  acid  transcripts.  Bioinformatics  is  then  used  to  extract  genomic  information. 

For  lithographically-fabricated  probe  arrays  (e.g.,  Affymetrix),  a  couple  of  VLSI-like  prob¬ 
lems  in  probe  placement  and  embedding  [5]  [6]  for  DNA  microarray  BioChips  design  have  been 
introduced.  However,  these  problems  are  often  dwarfed  by  the  complexity  of  the  more  general  bio¬ 
informatics  problem  of  probe  selection  (deciding  what  probes  should  be  placed  on  the  chip)  and 
image  processing  problem  of  the  chip  readout.  An  example  of  the  importance  of  image  processing 
can  be  garnered  from  Fig.  1(c),  which  shows  a  chip  image  after  the  unknown  bioagents  have 
reacted  with  the  probe  array.  An  alternative  to  lithographically  constructing  DNA  probes  directly  on 
the  chip  involves  synthesis  of  the  DNA  probes  off-chip,  and  then  inkjet  based  printing  [7]  or 
micro  fluidics-based  probe  transport  to  the  desired  location  on  the  chip  [8],  Micro  fluidics  is  also 
needed  for  BioChip  fluid  handling  in  high-throughput  screening  applications  [9], 

Due  to  the  combination  of  the  importance  of  microfluidics  platforms  and  the  similarity  of 
microfluidic  ICs  to  transistor  ICs  this  report  focuses  solely  on  microfluidics-based  Lab-on-a-Chips. 
LoC  design  issues  range  from  synthesis  and  physical  design  to  modeling  and  simulation-based  ver¬ 
ification.  The  next  section  in  this  report  describes  the  technology  and  design  issues  for  microfluid¬ 
ics-based  chips.  As  with  IC  design,  complexity  is  the  primary  motivating  factor  for  CAD,  as  will  be 
discussed  in  Section  2.2.  Top-down  design  methodologies  for  handling  this  complexity  are 
described  in  Section  3.1. 

2,1.  Lab-on-a-Chip 

pTAS  and  LoC  are  used  interchangeably  to  describe  chips  that  integrate  common  chemical 
or  biological  workbench  processes,  such  as  sample  pretreatment  and  transportation,  mixing,  reac¬ 
tion,  separation  and  detection  [10],  The  primary  advantage  of  miniaturization  advantages  include 
the  over  lOOOx  reduction  of  sample  consumption,  leading  to  improved  sensitivity  and  resolution  of 
biomolecular  detection.  Two  types  of  LoC  platfonns  are  being  developed,  categorized  by  their  fluid 
transport  mechanism. 
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2.1.1.  Digital  Microfluidics 

Transport  of  immiscible  gas  bubbles,  liquid  droplets  or  particles  (e.g.,  water  droplets 
immersed  in  oil  or  air)  have  been  explored  as  platforms  for  digital  micro  fluidics.  These  platforms 
can  be  used  as  a  discrete,  randomly  accessed,  multi-analyte  analyzer,  in  which  reagents  or  samples, 
compartmentalized  into  droplets,  can  be  processed  in  any  desired  order  in  parallel  fashion.  Of  the 
wide  variety  of  droplet/particle  transport  mechanisms  (usually  involving  surface  tension  modula¬ 
tion  [11]),  the  two  that  have  found  broad  biofluidic  applications  include  electrowetting  on  dielectric 
(EWOD)  and  dielectrophoresis  (DEP). 

In  EWOD,  the  electrolyte  droplet  is  placed  on  a  hydrophobic  surface.  This  surface  is  sepa¬ 
rated  from  the  electrodes  by  a  dielectric  layer,  as  shown  in  Fig.  2(a)  [12].  The  droplet  is  moved  by 
modulating  the  local  surface  tension  around  the  droplet  using  an  electric  field.  Droplets  can  be  cre¬ 
ated  from  a  reservoir  (microfluidic  digitization  —  Fig.  2(b)).  Additional  bio-processing  functions 
such  as  droplet  splitting  (Fig.  2(b)),  transportation  (Fig.  2(c)),  and  merging  (Fig.  2(d))  have  been 
successfully  demonstrated  [  1 3] [  14].  Recent  applications  include  calorimetric  enzyme -kinetic  assays 


Hydrophobic  layer 


Ground  electrode 


Hydrophobic 
dielectric  layer 

(a) 


^  Top  glass  plate 
Bottom  glass  plate 


Control  electrodes 


(b) 

/  '  /** — -i 

'  .  •*€  3* 

.  r  *) 

(c)  J 
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*  4 

(d) 
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FIGURE  2.  (a)  Cross  section  of  EWOD  device  [12];  (b)  Droplet  injection  from  reservoir  (left)  and  split 
(right);  (c)  Independent  and  simultaneous  droplet  transport;  and  (d)  two  droplets  are  merged  (left  and 
center),  and  then  merged  with  reservoir  (right)  1 14]. 
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[15]  and  on-chip  sample  processing  for  MALDI-Mass  Spectroscopy  [16]. 

DEP  forces  are  induced  on  polarizable  electrically  neutral  particles  due  to  an  induced  dipole 
when  subjected  to  a  non-unifonn  electric  field  [17].  Such  forces  have  been  used  widely  to  manipu¬ 
late  and  separate  cells  [18]  [19]. 

2.1.2.  Channel-based  Microfluidics 

Channel-based  micro  fluidics  [10]  predates  the  earliest  microarray  or  digital  micro  fluidics 
concepts,  and  is  still  more  widely  used  by  researchers  today  because  of  its  flexibility.  For  example, 
droplets  have  been  fonned  using  unstable  two  phase  flows  to  enhance  reaction.  In  terms  of  analo¬ 
gies  to  the  IC  world,  digital  microfluidics  tends  to  use  microarchitecture  abstractions,  while  chan¬ 
nel-based  microfluidics  tends  to  use  analog  circuit  abstractions. 

Channel-based  LoC  applications  range  from  analysis  of  amino  acids  [20]  and  DNA  [3] 
(exploiting  biochemical  separation),  to  chemical  synthesis  [21]  (exploiting  reaction).  Chips  inte¬ 
grating  both  reaction  and  separation  have  also  been  demonstrated  [22]. 

Fluid  flows  in  microchannels  exploit  two  primary  driving  forces:  pressure  and  electrokinetic 
(EK)  forces.  We  focus  on  EK  forces  here  because  (1)  they  can  be  created  by  applying  a  voltage  dif¬ 
ference,  like  transistor-based  ICs;  (2)  they  can  move  biomolecules  in  a  static  bulk  fluid,  needed  for 
molecular  separation;  and  (3)  they  tend  to  be  the  predominant  mechanism  for  complex  LoC 
designs. 

2.2.  LoC  Complexity 

Multiplex  complexity  arises  from  arrayed  integration  of  the  same  functional  subsystem.  If 


FIGURE  3.  (a)  Spiral  Capillary  Electrophoresis  separation  assay  for  amino  acids  |20];  (b)  Microreactors  for 
drug  discovery  [21];  (c)  Multifunction  immunoassay  that  integrates  mixing,  reaction,  and  separation  [22]. 
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the  unit  being  arrayed  is  identical,  e.g.,  the  chip  in  Fig.  1(b),  it  is  a  homogeneous  multiplex  (analo¬ 
gous  to  RAM  cells  being  arrayed  into  a  memory  IC).  If  the  units  have  different  sizes  but  the  same 
function,  then  it  is  a  heterogeneous  multiplex.  If  it  involves  different  functions,  such  as  when  mix¬ 
ing,  reaction  and  separation  are  integrated  as  in  the  immunoassay  in  Fig.  3(c),  then  it  is  a  multifunc¬ 
tion  system  (analogous  to  custom  ICs). 

Computational  fluid  dynamics  has  been  applied  to  channel-based  [23]  [24]  and  to  droplet- 
based  [25]  LoC  simulation  since  the  late-90’s.  This  detailed  simulation  provides  insight  at  the 
device  level  (e.g.,  for  a  single  turn  [26]  —  similar  to  device  simulation  in  transistor-based  ICs),  but 
is  impractical  for  system-level  simulation.  Typical  system-level  simulation  runs  can  take  days  and 
require  several  GB  of  RAM.  Additional  challenges  are  enumerated  in  [27].  While  fast  solvers  for 
gas  flow  have  been  developed  [28],  general  purpose  fast  solvers  are  still  not  available. 

3.  Methods,  Assumptions,  Procedures 

3.1.  Top-Down  Design  Methodologies 

To  overcome  these  complexity  limitations,  researchers  have  been  developing  top-down 
approaches  to  handle  complexity.  Hierarchy  is  used  to  overcome  the  limitations  of  continuum  simu¬ 
lations  for  channel-based  LoCs.  For  example,  the  immunoassay  in  Fig.  3(c)  can  be  partitioned  into 
mixing,  reaction,  injection  and  separation  subsystems.  We  call  this  the  functional  decomposition 
step  (breaking  the  system  into  its  functional  elements  similar  to  processor,  memory  and  I/O  in  com¬ 
puter  architecture).  Each  of  these  functional  elements  can  then  be  hierarchically  decomposed  into 
straight  channels,  elbows,  U-tums,  as  shown  by  the  schematic  in  Fig.  4.  This  is  the  geometrical 
decomposition  step  [29]. 

The  function-based  categorization  identifies  the  primary  physical  principles  that  need  to  be 
modeled  to  accurately  and  efficiently  capture  the  element  behavior.  The  geometry-based  categoriza¬ 
tion  allows  use  of  parametric  spatial  integration  to  obtain  geometrically  parameterizable  models, 
enabling  a  high  degree  of  model  reuse.  This  is  in  contrast  to  reduced-order  models  (ROMs)  which 
can  be  obtained  by  applying  nonlinear  circuit  model  reduction  techniques  [30].  These  ROMs  have 
the  potential  to  speed  up  simulation  of  a  specific  device,  but  need  to  be  regenerated  whenever  the 
geometry  changes,  which  is  undesirable  for  iterative  design. 

When  the  fluidic  network  only  involves  electroosmosis,  an  analogy  between  bulk  fluid  and 
current  flow  can  be  used  to  obtain  resistor-based  circuit  networks  for  fast  and  accurate  simulation 
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FIGURE  4.  System-level  Immunoassay  schematic  showing  straight  channels,  U  turns,  elbow  turns,  and  other 
elements  used  to  integrate  mixing,  reaction,  injection  and  separation  functions  on  a  single  LoC  [29]. 

[31].  This  fluid  transport  model  has  been  embellished  with  lumped  parameter  models  for  reaction, 
mixing,  and  separation  for  integrated  LoC  simulation  [32]. 

A  library  for  simulating  both  electrokinetic  and  pressure  driven  flows  in  microchannels  as 
well  as  droplet  transport  is  available  commercially  from  Coventor  as  a  fluidics  library  for  their 
ARCHITECT  simulator  [33]. 

These  first  generation  lumped  approaches  do  not  adequately  capture  dispersion  sources  in 
the  lumped  models.  Separate  continuum  simulations  are  needed  to  determine  the  correction  factors 
due  to  dispersion.  When  the  primary  dispersion  mechanism  is  geometrical,  this  implies  a  numerical 
simulation  run  for  every  change  in  design  geometry  [34],  making  these  approaches  untenable  for 
iterative  geometrical  design  for  LoC. 

To  overcome  these  limitations,  an  element  library  with  nine  mixing  elements,  ten  separation 
elements,  three  reaction  elements,  and  four  injection  elements  has  been  developed  at  Carnegie  Mel- 
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Ion  University.  This  library  has  been  used  to  create  schematic  networks  that  describe  dozens  of 
electrokine  tic  ally  driven  microfluidic  mixers,  separators,  and  integrated  assays  found  in  the  litera¬ 
ture  [35]. 


3.2.  Method  of  Moments 

Two  different  methods  were  used  to  arrive  at  closed-form  parameterized  models  for  the 
physiochemical  relationships  that  govern  behavior  in  biofluidic  microsystems.  In  essence  these 
approaches  aimed  at  replacing  the  partial  differential  governing  equation.  The  first  of  these  methods 
is  the  Method  of  Moments  and  is  described  here. 

The  species  band  concentration  c(y,z,t )  within  a  separation  channel  is  governed  by  the  con¬ 
vection-diffusion  equation  [36] 


dc  dc  ^ 

dt  dz 


f  d2c  d2c^ 
dz 2  dy2 


(1) 

y  \jy  j 

where  y  and  z  are  the  transverse  and  axial  positions  respectively,  t  is  the  separation  time  from  the 
channel  entrance.  The  width  of  the  species  band  can  be  characterized  by  variance,  the  square  of  the 
standard  deviation  of  the  average  concentration  profile  cm,  which  is  defined  as 


a2  = 


f  ( z~z)2cm-dz 

J  —  CO 


•  dz 


(2) 


where  y  is  the  axial  position  of  the  species  band’s  centroid  in  the  channel. 

Eq.  (1)  can  be  reformulated  into  a  more  tractable,  reduced-dimension  form  in  terms  of  spa¬ 
tial  moments  of  the  species  concentration.  Such  moments  are  capable  of  describing  the  species 
band’s  main  characteristics  such  as  mass  distribution,  skew  and  variance  without  solving  for 
detailed  concentration  distributions.  We  introduce  a  new  coordinate  frame,  which  moves  at  the  spe¬ 
cies  band’s  average  velocity  U,  and  normalize  the  equation  to  reduce  all  variables  into  dimension¬ 
less  forms.  Define  a  dimensionless  axial  coordinate  g,  transverse  coordinate  r/  and  time  rby 


4  =  {z  —  Ut)/w,  T]=y/w,  r  =  Dt/w 2  (3) 

In  tenns  of  these  dimensionless  variables,  Eq.  (1)  is  rephrased  to  the  following  form  in  the 

concentration  c{^,ij,r) : 


dc  d2c  d2c  _  dc 


(4) 
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where  P e  =  Uw/D  is  the  Peclet  number  representing  the  ratio  of  convection  and  diffusive  transport 
rates,  and  x  is  the  normalized  species  velocity  relative  to  the  average: 

j(?7)  =  (m -£/)/£/  (5) 

We  now  recast  Eq.  (4)  in  tenns  of  spatial  moments  of  the  species  concentration.  If  the  spe¬ 
cies  band  is  entirely  contained  in  the  channel,  Eq.  (4)  holds  effectively  over  the  axial  domain 
— oo  <  S,  <  oo  (the  transverse  domain  is  0  <  77  <  1 ),  such  that  c— » 0  as  ^—>±00.  Therefore,  we  can  define 
spatial  moments  of  the  species  concentration  by 

cP(Jl,*)  =  \"jPc{Z,il,'c)dZ,  mp  (r)  =  f'0cpd/)  (6) 

Here,  cp  is  the  /fh  moment  of  the  species  concentration  in  the  axial  filament  at  //,  and  mp  is 

the  /ith  moment  of  the  average  concentration  of  the  band,  respectively.  Note  that  as  a  consequence  of 
the  coordinate  transformation  Eq.  (3),  all  moments  are  taken  with  respect  to  the  moving  frame 
(£  77).  For  purposes  of  simulating  analyte  dispersion,  it  suffices  to  obtain  the  moments  up  to  the  sec¬ 
ond  order.  Specifically,  cq  provides  the  transverse  distribution  of  the  species  mass  in  each  axial  fila¬ 
ment  within  the  channel  and  m0  is  the  total  species  mass  and  can  be  chosen  as  m()=  1  without  losing 
generality.  Next,  cj  gives  the  axial  locations  of  the  centroid  of  the  axial  filaments  in  the  species 
band,  and  hence  indicates  the  skew  of  the  band.  Then,  ni\,  the  widthwise  average  of  cj,  is  the  axial 
location  of  the  centroid  of  the  entire  species  band  in  the  frame  (£  rj)  and  always  zero  for  this  study 
[37].  Finally,  m2  can  be  used  to  detennine  the  variance  cr  of  the  species  band  by 

a1  =  w2  (m2  / m0  -  m2  / t/Zq  ) 

3.3.  Reduced  Order  Modeling 

The  second  approach  to  getting  parameterized  models  is  reduced  order  modeling.  The  tree  in 
Fig.  5  distinguishes  between  numerical  and  analytical  modeling  methods.  When  problems  cannot 
be  solved  analytically  with  closed-fonn  solutions  to  the  governing  partial  differential  equations, 
additional  approximations  must  be  made.  Numerical  modeling  makes  these  approximations  in  vari¬ 
ous  ways.  In  full  order  methods,  the  original  equations  are  discretized  and  solved  in  full  form.  This 
new  system  of  equations  can  be  very  large  and  extremely  computationally  expensive  to  evaluate. 
Alternate  numerical  methods  seek  to  reduce  this  computational  complexity  by  removing  nonessen¬ 
tial  degrees  of  freedom.  These  so  called  ‘reduced-order’  methods,  such  as  state-space  reductions, 
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FIGURE  5.  Related  modeling  methods. 

reduce  the  scale  of  the  discretization  and  eliminate  unnecessary  degrees  of  freedom  in  the  discreti¬ 
zation  using  basis  reduction  methods  [38].  In  other  reduced  order  methods,  if  specific  system 
responses  are  desired,  such  as  the  yield  of  a  chemical  reaction,  a  single  function  can  be  found 
through  parameter  estimation  of  a  known  functional  form  obtained  from  physical  knowledge  of  the 
system  [39].  Alternatively,  if  the  system  is  a  ‘black-box’,  the  basis  of  its  functional  representation  is 
unknown  due  to  a  lack  of  physical  knowledge.  In  these  cases,  non-linear  regression  techniques, 
such  as  the  neural  network  can  be  used  where  the  choice  of  functional  form  is  not  necessary 
[40]  [41]  [42].  The  focus  of  the  modeling  in  this  work  is  to  create  response  surface  models  targeted  at 
design  synthesis.  The  developed  models  are  valid  only  within  the  specified  design  space,  which 
may  or  may  not  be  rectangular. 

3.3.1.  Neural  Network  Reduced  Order  Model 

The  creation  of  a  model  using  numerical  simulations  involves  the  definition  of  a  mapping 

between  the  model  inputs  and  parameterized  outputs.  The  numerical  simulation,  which  approxi¬ 
mates  the  physics  of  the  real  model,  is  used  to  sample  and  define  this  mapping,  followed  by  a  multi¬ 
dimensional  regression  to  create  a  closed-form  model.  This  model  approximates  the  device's  actual 
performance  mapping  function  (the  real  model).  A  neural  network  is  a  robust  regression  tool  for 
identifying  this  mapping  because  it  is  efficient  with  sparsely  sampled  data,  can  identify  mappings 
without  any  prior  information  about  the  desired  fonn  of  the  model  or  physics,  and  results  in  a 
closed-form  function  parameterized  by  the  weights  inherent  to  the  connections  between  the  nodes 
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of  the  NN.  In  this  work  all  calculations  are  performed  in  MATLAB  [43]. 

When  the  nodes  of  a  NN  are  arranged  in  layers,  with  connections  only  between  adjacent  lay¬ 
ers,  the  topology  is  referred  to  as  feed-forward.  An  example  of  this  topology  is  shown  in  Fig.  6.  The 
NN  operates  by  passing  the  data  at  its  input  nodes  through  a  series  of  weights  to  the  nodes  in  the 
layers  above.  Each  node  in  the  layer  takes  the  weighted  sum  of  the  output  of  the  nodes  below  as  the 
argument  of  some  arbitrary  bound  non-linear,  or  unbound  linear  transfer  function.  The  use  of  non¬ 
linear  transfer  functions  in  all  nodes  between  the  input  and  output  layers  and  linear  transfer  func¬ 
tions  at  the  output  nodes  gives  the  NN  the  useful  property  of  being  a  universal  function  approxima¬ 
tor  when  at  least  two  of  these  hidden  layers  are  present,  and  sometimes  only  one  layer  when  the 
approximated  function  is  bound  [44].  This  means  given  enough  nodes  in  the  hidden  layers,  the  NN 
will  be  able  to  approximate  a  function  of  any  degree  of  non-linearity,  with  a  finite  number  of  dis¬ 
continuities,  whether  bound  or  otherwise,  to  an  arbitrary  degree  of  accuracy  [45].  So  with  a  certain 
degree  of  confidence,  it  can  be  assumed  that  the  NN  regressor  will  provide  a  robust  method  for 
mapping  the  perfonnance  function  for  a  very  large  class  of  device  models,  given  sufficient  sam¬ 
pling  data  to  represent  the  details  of  the  underlying  structure  of  the  mapping. 

To  uniquely  define  the  NN  model,  the  number  of  layers,  nodes  per  layer,  node  transfer  func¬ 
tion,  and  interconnecting  weights  must  be  specified.  In  general,  the  number  of  layers  and  node 
transfer  functions  are  chosen  to  achieve  the  universal  approximator  criterion,  namely  two  hidden 
layers  with  bound  non-linear  transfer  functions  and  linear  transfer  functions  on  the  output  nodes.  In 


FIGURE  6.  (a)  The  activation,  of  a  node  is  a  function  of  the  weighted  inputs  of  the  previous  layer's 
activations,  (b)  General  feedforward  neural  network  topology  utilizing  multiple  hidden  layers  with  an 
unspecified  number  of  nodes  with  non-linear  activation  functions,  and  a  single  output  node  using  a  linear 
activation  function. 
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some  of  the  models,  namely  for  the  injectors,  only  a  single  hidden  layer  of  nodes  is  required  since 
the  targeted  performance  functions  are  almost  strictly  monotonic  and  bound.  The  weights  connect¬ 
ing  the  nodes  are  chosen  so  as  to  minimize  an  error  criterion,  such  as  the  mean  squared  error  of  the 
NN  model  to  the  sampled  mapping  data.  This  process  is  called  training  the  neural  network.  To 
choose  the  number  of  nodes  in  the  hidden  layers,  the  number  of  nodes  is  varied  while  the  NN 
trained.  The  number  of  nodes  which  minimizes  the  degree  of  underfilling  and  overfitting  is  chosen 
for  the  final  model. 

Underfitting  and  overfitting  manifests  itself  through  the  bias  and  variance  of  the  network  as 
a  function  of  the  complexity  of  the  network  (more  nodes  makes  a  more  complicated  network).  The 
mean  squared  error  decomposition  for  a  particular  data  point,  x,  for  noiseless  source  data  leads  to 
the  following  metric  for  the  networks  estimate: 

MSE(f(x ))  =  E[(f(x)-E\f(x)])2]  +  (E\f(x)]-f(x))2  (7) 

where  f(x)  is  the  network's  estimate  for  the  mapping  at  point  x,  f(x)  is  the  actual  target  data  at  that 

point,  and  E[...\  represents  the  expectation  for  an  ensemble  of  training  data  sets.  The  first  summand 
is  a  measure  of  the  variance  of  the  estimate  at  that  point,  and  the  second  summand  is  a  squared  mea¬ 
sure  of  the  bias,  or  difference  between  the  expected  estimate  and  the  actual  data  at  that  point,  for  sets 
of  training  data  [46]  [47]. 

Networks  of  lower  complexity  (fewer  nodes)  will  tend  to  have  a  high  bias  and  a  low  vari¬ 
ance,  whereas  networks  with  high  complexity  (more  nodes)  will  tend  to  have  low  bias  and  high 
variance.  To  find  the  optimal  complexity  for  a  model,  the  number  of  nodes  in  the  hidden  layer  is 
chosen  so  as  to  find  the  best  trade-off  between  bias  and  variance.  Two  types  of  error  measurements 
are  conducted  during  the  construction  of  a  NN:  l)Training  error  measures  the  error  of  the  network 
to  the  original  data  set  (the  training  data  set).  Since  the  training  process  is  an  optimization  problem 
that  constrains  the  network  to  have  a  low  bias  on  this  data  set,  this  is  not  a  good  indicator  of  the 
overall  generalized  error  of  the  network  (for  example,  it  doesn't  indicate  the  bias  or  variability  of  the 
network  to  data  it  has  not  seen  in  the  training  data  set).  2)To  measure  the  variance  of  the  network,  a 
different  test  data  set  must  be  used,  consisting  of  points  not  included  in  the  training  data  set,  which 
provides  a  measure  of  the  variability  at  points  where  the  network  was  not  constrained  during  its 
learning  process.  This  is  a  more  generalized  error  as  indicated  by  Eq.  (7).  A  network  with  a  high 
bias  or  variance  will  produce  a  large  generalization  error,  and  this  is  the  error  that  must  be  mini- 
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mized  to  choose  the  optimal  complexity  of  a  network.  Efficient  methods  for  measuring  this  general¬ 
ization  error  on  test  data  sets  is  outline  in  the  next  section. 


3.3.2.  Model  Selection 

.A  technique  referred  to  as  k-fold  cross  validation  (KFCV)  is  a  practical  method  for  measur¬ 
ing  the  generalized  error  in  Eq.  (7).  This  technique  is  efficient  since  it  does  not  require  a  separate 
validation  set.  All  data  is  used  in  the  training  of  the  neural  network.  As  shown  in  Fig.  7,  this  method 
randomly  sorts  the  training  data  set  into  k  subsets  (folds).  The  network  is  trained  on  a  set  consisting 
of  k- 1  subsets  and  trained  on  the  one  remaining  subset.  In  order  to  avoid  biasing,  each  subset  plays 
the  role  of  the  validation  set  exactly  once.  The  errors  measured  for  each  validation  subset  iteration 
are  then  averaged  together  to  provide  an  unbiased  measure  of  the  generalized  error  It  has  been 
detennined  that  k  of  5  or  10  provides  the  best  estimate  of  the  generalization  error  [38][39]. 

Other  types  of  validation  include,  the  simple  test  set  validation.  In  this  method,  a  fixed  num¬ 
ber  of  data  points  are  removed  from  the  acquired  training  data  and  are  used  after  training  to  measure 
the  ability  of  the  model  to  generalize  to  unknown  data.  This  method  is  undesirable  for  several  rea¬ 
sons.  First,  it  reduces  the  amount  of  data  available  for  training,  which  is  important  when  the  data  set 
is  not  large  to  start.  Second,  it  has  a  large  bias  error  since  the  removed  data  is  a  fixed  set  that  clearly 
only  significantly  affects  the  accuracy  of  the  model  in  the  region  of  missing  data.  It  is  also  unclear 
how  to  decide  which  data  points  should  be  removed  from  the  original  set. 

Leave-one-out  cross  validation  improves  on  the  test  set  validation  method  by  training  a 


random  sort  S 
for  i  =  1 : k 
D±  =  S (1 :N/k) 
end 

for  i  =  1 : k 
T  =  { } 

V  =  {  D±  } 

NN  =  CTrain  NN  on  T> 
err±  =  MSE (NN (V) ) 
end 


S :  N  points 

{ D  |  :  k  subsets 

k- 1  training  folds 
1  validation  fold 


□-o 

enq 


err  =  mean({erri}) 


□ 


FIGURE  7.  k-fold  cross  validation  algorithm.  First  divide  the  N  points  in  the  dataset  into  k  subsets  (3  to  5 
subsets  is  typical).  Then,  identify  one  subset  (  i  )  as  the  validation  set,  used  to  compare  with  the  trained 
neural  network.  The  other  k- 1  subsets  ( ~i ),  are  used  to  train  the  neural  network.  Each  subset  is  the 
validation  set  exactly  once.  The  error  for  each  iteration  is  averaged  at  the  end  to  provide  an  estimate  of  the 
generalization  error  for  the  network.  This  is  an  efficient  use  of  data,  since  all  data  is  used  for  training  and 
validation.  There  is  no  exclusive  validation  set. 
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model  N  times  for  a  data  set  with  N  points,  leaving  one  point  out  of  the  training  data  set  each  time, 
and  measuring  the  error  of  the  model  on  that  point  afterwards,  [41].  This  removes  the  bias  intro¬ 
duced  from  using  a  fixed  test  set,  but  since  each  of  the  N  training  data  sets  are  very  similar  since 
only  one  point  is  missing,  the  variance  of  the  resulting  model  can  be  very  high. 


3.3.3.  Design  of  Experiments 

The  neural  network  can  only  be  as  accurate  as  the  data  used  for  training.  There  are  several 
factors  to  consider  when  choosing  data  to  train  the  network,  these  and  related  issues  are  generally 
described  as  design  of  experiments  or  data  mining  [46].  First  is  the  number  of  total  data  points  to 
provide.  There  is  no  fixed  number  of  samples  that  can  be  used  in  general,  since  it  depends  on  the 
number  of  dimensions  and  the  complexity  of  the  function  in  the  space.  However,  rapidly  changing 
functions  in  high  dimensions  will  require  many  more  samples  than  slowly  changing  functions  in 
low  dimensions.  As  an  example,  for  the  injector  models  in  the  following  section  with  four  dimen¬ 
sions,  on  the  order  of  1000  points  are  used.  Second  is  where  to  sample  in  the  space  spanning  the 
dimensions  of  the  function  to  be  learned.  If  knowledge  of  the  underlying  function  shape  is  used, 
samples  can  be  increased  near  points  of  rapid  change.  Otherwise,  uniform  sampling  throughout  the 
space  is  used.  In  order  to  maximize  the  number  of  samples  per  dimension,  it  is  best  to  sample  the 
data  not  on  a  unifonn  grid.  Fig.  8  shows  an  example  in  two  dimensions  where  the  number  of  sam¬ 


ples,  ns,  is  taken  on-  and  off-  grid.  When  taken  on- grid  there  are  fns  unique  points  per  dimension 
(in  an  k-dimensional  space),  where  off-grid  results  in  ns  unique  points  per  dimension. 


The  off-grid  results  can  be  generated  using  pseudo-random  sequence  generators,  which 
asymptotically  have  a  lower  discrepancy  (defined  in  a  moment)  than  purely  random  sequences.  The 
benefit  of  using  a  grid  to  sample  is  that  the  points  are  equidistantly  separated  eliminating  all  bias  in 
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FIGURE  8.  (a)  Uniform  sampling  on  a  grid  in  two-dimensions  with  a  total  of  16  data  points  resulting  in 
effectively  4  data  points  per  dimension,  (b)  Uniform  sampling  off-grid  in  two  dimensions,  again  with  16 
data  points,  resulting  in  effectively  16  points  per  dimension,  a  better  more  denser  sampling  per  dimension. 
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the  sampling  sequence  —  one  area  isn’t  over-sampled  while  another  area  with  important  structure 
might  be  undersampled.  Unfortunately,  as  mentioned  before,  the  use  of  a  grid  reduces  the  number 
of  unique  samples  per  dimension.  This  trade-off  can  be  balanced  by  using  a  low-discrepancy 
sequence  that  maximizes  the  sampling  per  dimension,  while  minimizing  the  clustering  of  high-dis¬ 
crepancy  sequences.  These  low-discrepancy  sequences  can  be  achieved  asymptotically  using 
pseudo-random  sequences,  such  as  the  common  Sobol  [48][49],  Faurer  [50],  Halton  [51],  Hainmer- 
sley  [52],  or  Niederreiter  [53]  generators.  The  quasi-random  sequences  were  originally  created  to 
make  more  rapidly  converging  high-dimensional  integrals  using  Monte-Carlo  evaluation,  but  the 
results  apply  to  neural  networks,  as  well. 

Discrepancy  can  be  defined  on  the  unit  hypercube  as  the  maximum  difference  between  the 
volume  (for  any  volume  defined  within  the  hypercube  —  with  one  vertex  at  the  origin  for  the  star- 
discrepancy)  and  the  average  number  of  points  included  within  the  volume, 


J  V 

D  =  sup 

M  e  [0,\)d 


i  N  d  d 

mZ  II"  II' 

i=y=  1  7=1 


(8) 


where  0  is  the  unit  step  function,  n/7  is  the  ith  sampled  point  out  of  N  in  the  jth  dimension  out  of  d, 
and  Xj  represents  the  edge  in  the  jt!l  dimension  of  the  variable  cube  measuring  enclosed  points  within 
the  unit  hypercube.  A  sequence  {a  A  is  considered  uniformly  distributed  if  the  average  number  of 
points  within  a  volume  equals  the  size  of  the  volume  as  the  sequence  becomes  infinite  (i.e.  uniform 
distributions  have  all  possible  volumes  “completely”  full  with  points,  whereas  non-unifonn  distribu¬ 
tions  will  leave  some  volumes  more,  or  less,  full). 


For  this  work,  the  non-asymptotic  nature  of  the  datasets,  with  on  the  order  of  1000  points, 
means  the  discrepancy  of  pseudo-random  sequences  are  not  significantly  better  than  a  randomly 
perturbed  deterministic  uniform  sequence.  Further,  the  randomly  perturbed  deterministic  sequence 
has  several  additional  post-processing  benefits  not  possible  with  quasi-random  sequence.  For  exam¬ 
ple,  if  certain  regions  of  the  problem  domain  are  marked  as  infeasible,  the  perturbed  sequence  can 
be  generated  so  as  to  “push”  points  out  of  the  infeasible  region.  This  prevents  running  simulations 
on  points  that  will  eventually  be  excluded  from  the  model.  Such  perturbations  can  be  applied  to  a 
quasi-random  sequence,  but  only  after  the  sequence  has  been  generated  as  a  whole,  thus  such  per¬ 
turbations  negate  the  low-discrepancy  benefits  of  using  such  sequences  in  the  first  place. 

The  improvement  in  training  accuracy  for  scattered  data  points  rather  than  uniform  gridded 
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FIGURE  9.  A  neural  network  is  trained  on  the  function  shown  on  the  left.  The  network  has  a  single  hidden 
layer  of  tansig  neurons  and  one  linear  output  layer.  It  is  trained  for  450  max  iterations,  and  the  sample  size 
is  varied.  The  error  measured  is  the  mean  squared  error  of  the  network  to  the  training  data  (not  a  measure 
of  the  generalization  error).  The  Sobol  quasirandom  sequence  is  the  best  performing,  followed  by  the 
Gaussian  sequence,  and  finally  the  uniform  grid.  The  Gaussian  sequence  is  a  simpler  algorithm  and  is 
compatible  with  domain  shaping  methods. 


points,  can  be  shown  with  an  example.  In  this  example,  three  methods  are  compared  for  training  a 
NN  on  a  two-dimensional  function: 


f(x,y)  =  cos(x)sin(v)exp 


-x-y 
2n 


(9) 


First,  x  and  v  are  sampled  from  0  to  2 n  on  a  uniform  grid,  using  the  meshgridQ  command  in 
MATLAB.  Then  x  and  y  are  sampled  from  0  to  2n  with  each  point  perturbed  off  the  uniform  grid  by 
an  amount  determined  by  a  gaussian  with  a  standard  deviation  equal  to  3%  of  the  domain  size  (27t). 
In  this  case,  if  any  points  go  outside  the  defined  domain,  they  are  resampled  until  they  fall  within  0 
to  2n.  Finally,  x  and  v  are  sampled  from  0  to  2n  using  the  Sobol  pseudo  random  sequence  generator 
[48] [49]  [54],  obtained  from  [55].  Using  each  of  these  sampling  methods,  a  NN  is  trained  and  the 
error  to  the  training  data  is  measured  and  plotted  in  Fig.  9.  The  neural  network  has  the  topology  of 
Fig.  6b  with  a  single  hidden  layer  of  10  neurons  with  a  tansig  transfer  function, 


fix) 


1 +  exp(-2x) 


-  1 


(10) 


19 


and  a  single  output  layer  with  a  linear  transfer  function, 

fix)  =  x  (11) 

The  network  is  trained  for  a  maximum  of  450  iterations  using  the  Levenberg-Marquardt 

algorithm  [56].  Since  the  initialization  of  the  algorithm  is  random,  each  point  in  the  graph  of  Fig.  9 
represents  the  average  of  40  repeated  training  experiments.  The  initial  slope  in  the  individual  traces 
is  most  likely  due  to  the  oversimplified  representation  of  Eq.  (9)  from  undersampling  with  the 
smaller  datasets.  As  the  dataset  increases,  the  additional  structure  of  Eq.  (9)  makes  the  training 
more  difficult  leading  to  the  initial  rise  in  error  as  it  progresses  to  some  other  asymptotic  value  that 
depends  on  the  complexity  of  the  network.  It  can  be  concluded  from  Fig.  9  that  the  Sobol  psuedo 
random  sequence  provides  the  best  results,  followed  by  the  Gaussian  sequence,  and  finally  the  uni¬ 
form  sequence.  Due  to  the  more  complex  algorithm  and  loss  of  the  low-discrepancy  qualities  when 
the  points  are  moved  to  avoid  infeasible  areas  of  the  parameter  space,  the  Gaussian  algorithm  will 
be  the  default  algorithm  used,  unless  otherwise  noted. 

4.  Results,  Discussion 

4.1.  Component  Models  Development  and  Verification 

The  goal  of  each  behavioral  model  is  to  capture  the  input-output  signal  flow  relationship  of 

the  pin  values  that  define  biofluidic  state  at  the  inlet  and  outlet  of  each  element  [57].  This  captures 
the  physical  phenomena  being  modeled  in  that  element.  In  addition,  an  electrical  resistance  is  asso¬ 
ciated  with  each  element  to  relate  the  EK  current  flow  through  the  element  to  the  inlet  and  outlet 
voltages.  In  contrast  to  the  bottom-up  reduced-order  model  approaches,  our  behavioral  models  pos¬ 
sess  several  important  attributes  to  enable  accurate  and  efficient  system-level  simulations  of  com¬ 
plex  LoCs.  Our  analytical  models  effectively  account  for  the  same  multi-physics  (e.g., 
electrostatics,  fluidics  and  mass  transfer)  as  numerical  simulation  tools.  They  do  not  require  any 
parameters  from  user-conducted  experiments  or  numerical  simulations  to  capture  interactions 
between  the  elements,  and  hence  provide  seamless  model  interconnectivity.  Most  importantly,  they 
are  in  closed-form  and  are  all  parameterized  by  element  dimensions  and  material  properties;  there¬ 
fore  they  are  reusable,  fast  to  evaluate,  and  well  suited  for  an  iterative  simulation-based  design 
methodology. 

As  discussed  above,  depending  on  the  physical  phenomena  of  individual  devices,  contents 
of  the  behavioral  model  libraries  will  be  different.  Hence,  models  for  the  micromixers  and  electro- 
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phoretic  separation  systems  will  be  developed  separately,  and  are  available  in  separate  model  librar¬ 
ies  for  schematic-based  simulations. 


4.1.1.  Electro  kinetic  Passive  Micromixers 

The  EK  passive  micromixer  library  consists  of  models  for  nine  elements,  which  includes 
reservoirs  (sample  and  waste),  slightly  tapered  straight  mixing  channel,  turns  (90°  or  180°,  clock¬ 
wise  or  counter-clockwise),  as  well  as  converging  and  diverging  intersections.  In  this  section,  we 
will  present  behavioral  models  for  basic  elements  such  as  the  slightly  tapered  mixing  channel,  con¬ 
verging  and  diverging  intersections.  Other  elements  can  be  modeled  in  a  similar  fashion. 

Slightly  Tapered  Straight  Mixing  Channel.  The  tapered  straight  mixing  channel,  in  which  differ¬ 
ent  samples  a  and  b  mix  with  each  other,  has  one  inlet  and  one  outlet,  with  different  cross-sectional 
area.  It  is  critical  in  designing  a  geometrical  focusing  micromixer  [58],  Electrically,  it  is  modeled  as 
a  resistor  and  the  resistance  is  given  by 

D  O  dZ 

w(z)h(z)Ce  (12) 

where  w  and  h  are  the  channel  width  and  depth  (both  are  functions  of  the  axial  coordinate  z),  Ce  is 
the  electric  conductivity  of  the  buffer  solution  in  the  channel.  As  a  special  case,  within  a  straight 
channel  with  the  uniform  cross-section,  Eq.  (12)  can  be  reduced  to 


R  = 


L 

hwCe 


(13) 


To  obtain  the  sample  concentration  profile  at  the  outlet,  we  partition  the  slightly  tapered 
straight  channel  into  a  series  of  segments  (segment  number  tends  to  infinity),  each  with  unifonn 
cross  section.  In  each  segment,  the  convection-diffusion  equation  is  solved  to  establish  the  input- 
output  relationship  of  concentration  coefficients  between  the  inlet  and  outlet  of  the  segment.  Then 


all  the  segmental  solutions  are  multiplied  and  the  concentration  coefficients  (n  =  0,1,2... )  at  the 
channel  outlet  are  attained  as  [59] 


2  2  LD 
-yn  n 


d{oul)  =  d(in)e  E,„m 


(14) 


where  d(‘n) ,  wjn  and  Ein  are  the  concentration  coefficients,  channel  width  and  electric  field  at  the  inlet, 
y  is  a  factor  capturing  the  effect  of  the  cross-sectional  shape  on  mixing  [59],  D  and  //  are  the  diffu- 
sivity  and  electrokinetic  (including  both  electroosmotic  and  electrophoretic)  mobility  of  the  sample 
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FIGURE  10.  (a)  Behavioral  model  structure  for  the  converging  intersection  in  the  micro  mixer.  Index  /,  r  and 
out  represent  the  quantities  at  the  left  inlet,  right  inlet  and  outlet.  Both  electric  ( V)  and  biofluidic  (</,„)  pins 
are  defined  at  the  terminals  of  the  model.  Electrically,  it  is  modeled  as  a  combination  of  three  resistors  (R/,Rr 
and  Rout )  with  zero  resistance.  Different  sample  concentration  profiles,  c^rj)  and  cr(rj),  at  inlets  are  merged 
and  compressed  at  the  outlet  cout(if).  (b)  Behavioral  model  structure  for  the  diverging  intersection  in  the 
micromixer.  Similarly,  index  /,  r  and  in  represent  the  quantities  at  the  left  outlet,  right  outlet  and  inlet. 
Sample  concentration  profile  at  inlet,  c,„( q)  is  split  and  stretched  out  into  two  parts,  c^rj)  and  cr(r/),  flowing 
out  of  the  outlets. 

and  L  is  the  channel  length.  The  special  case  of  a  straight  channel  with  the  unifonn  cross-section 
yields  y  =  1. 

Converging  Intersections.  Fig.  10  shows  the  behavioral  model  structure  of  converging  and 
diverging  intersections  used  in  micromixers  [60].  Arrows  at  pins  indicate  the  direction  of  signal 
flow  for  computing  biofluidic  pin  values  and  state.  The  converging  intersection  has  two  inlets  and 
one  outlet,  and  acts  as  a  combiner  to  align  and  compress  upstream  sample  streams  of  an  arbitrary 
flow  ratio  s  (defined  below)  and  concentration  profiles  side-by-side  at  its  outlet  (Fig.  10a).  As  its 
flow  path  lengths  are  negligibly  small  compared  with  those  of  mixing  channels,  such  an  element 
can  be  assumed  to  have  zero  physical  size,  and  electrically  represented  as  three  resistors  with  zero 
resistances  between  each  terminal  and  the  internal  node  N, 

Rl=K=R0u,=  0  (15) 

Flere,  N  corresponds  to  the  intersection  of  flow  paths  and  the  subscripts  /,  r  and  out  represent 

the  left  and  right  inlets,  and  the  outlet,  respectively.  Defining  <$  and  d^'}  (m= 0,1,2...)  as  the  Fou¬ 
rier  coefficients  of  the  sample  concentration  profiles  at  the  left  and  right  inlets  respectively,  then  the 

coefficients  d[out^  («  =  0,1,2...)  of  the  profile  at  the  outlet  (cout  (77))  are  related  to  d$  and  d^  by, 
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S  <  1J  <  1 


Eq.  (16)  shows  that  the  concentration  profile  at  the  outlet  can  be  treated  as  a  superposition  of 
the  scaled-down  profiles  from  both  inlets,  where  s  =  ql/(qI+qr)  =  I,/(ll+Ir)  denotes  the  interface 

position  (or  flow  ratio,  the  ratio  of  the  left  flow  rate  <27  to  the  total  flow  rate  (q,  +  qr) )  between 
incoming  streams  in  the  normalized  coordinate  at  the  outlet  (note  that  flow  rates  q/  and  qr  are 
respectively  linear  with  the  electric  currents  //  and  I 

Solving  Eq.  (16)  yields  d[°“1'1  as, 


diout)=4l)s  +  4r\l-s) 


,  s  m^ns  f  _ •  /  r  \  ,  r  _•  /  r  \  m=ns  m-n(\  s ) 
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m= 0 


cos  (E2/2)  sin  (Fj/2)  cos(/7j/2)sin(F2/2) 

^  +  7^ 


where  fx=(m-ns)n,  f2=(m  +  ns)x,  Fl={m  +  n-ns)x  and  /f  =  +  .  Since  the  sample  concen¬ 

tration  profiles  at  the  inlets  are  scaled  down,  the  Fourier  series  mode  at  the  inlets  is  not  orthogonal 
to  that  at  the  outlet.  Therefore,  the  calculation  for  the  coefficient  of  a  certain  Fourier  mode  at  the 
outlet  also  depends  on  the  other  modes  at  the  inlets. 


Diverging  Intersections.  The  diverging  intersection  has  one  inlet  and  two  outlets  and  is  the  dual  of 
the  converging  intersection.  It  splits  the  incoming  flow  and  electric  current  into  two  streams  that 
exit  from  the  outlets.  It  can  also  be  represented  by  three  zero-resistance  resistors, 


R.n=Rl=Rr=0  (18) 

where  subscripts  in,  l  and  r  represent  quantities  at  the  inlet,  the  left  and  right  outlets. 

Defining  (m=0, 1,2 _ )  as  the  Fourier  coefficients  of  the  sample  concentration  profile  at 

the  inlet.  Then  the  coefficients  at  the  left  and  right  outlets  and  d\’ 1  arc  respectively  given  by 
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(19) 


and 


oo 
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-n/s 
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m#n/(l—s ) 

^l>o=2  Z  4wV2sin(^)M+  Z  (-i)"""4'") 

/w=0  m—0 

where  fl=(n-ms)K,  f2=[n  +  ms)7r,  Fx  =  {n  +  m-ms)K,  F2  =  {n-m  +  ms)K,  <j>x=msn  and  $  =  w(l-.s);r . 
Similar  to  the  converging  intersection,  5  is  the  normalized  splitting  position  (or  ratio). 


It  should  be  pointed  out  that  in  contrast  to  the  resistor-based  mixing  models  [32]  [60]  that 
take  advantage  of  the  analogy  between  fluidic  and  sample  transports  and  only  convey  the  average 
concentration  values  through  the  entire  network,  our  models  (Eqs.  (14),  (17),  (19)  and  (20))  propa¬ 
gate  sample  concentration  profiles  characterized  by  the  Fourier  series  coefficients.  This  removes  the 
requirement  of  complete  mixing  (along  channel  width)  at  the  end  of  each  channel  in  the  network 
imposed  by  the  resistor-based  model  and  allows  for  optimal  design  of  both  effective  and  efficient 
micromixers. 


4,1.2.  Electrophoretic  Separation  Chips 

The  electrophoretic  separation  library  includes  models  for  ten  basic  elements:  turns  (90°  or 
180°,  clockwise  or  counter-clockwise),  straight  channel,  detector,  injector,  injection  channel  and 
reservoirs  (sample  and  waste).  In  this  section,  behavioral  models  for  basic  elements  such  as  separa¬ 
tion  channels  (straight  and  turn)  will  be  developed  to  analyze  the  band-spreading  effect  caused  by 
molecular  diffusion  and  turn  dispersion.  Additionally,  a  detector  model  applicable  for  both  DC  and 
transient  analysis  will  be  presented.  Models  of  the  other  elements  can  be  derived  using  the  same 
principles. 

Fig.  11  shows  the  behavioral  model  structure  of  electrophoretic  separation  channels  (straight 
or  turn).  Arrows  indicate  the  direction  of  signal  flow  for  calculating  biofluidic  pin  values  and  state. 
Electrically,  separation  channels  are  modeled  as  resistors  in  the  same  way  as  the  unifonn  straight 
mixing  channels  (for  a  constant-radius  turn,  L  in  Eq.  (13)  is  replaced  by  L  =  rc0  ?  where  rc  and  6  are 
the  mean  radius  and  angle  included  by  the  turn,  see  [37]  [59]  for  detailed  geometrical  interpretation). 
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FIGURE  11.  Behavioral  model  structure  for  separation  channels  in  electrophoretic  separation  microchips. 
Index  in  and  out  represents  the  quantities  at  the  inlet  and  outlet.  Both  electric  (V)  and  biofluidic  ( t ,  Sn  ,  a2 
and  A)  pins  are  defined  at  the  terminals  of  the  model.  Electrically,  the  channel  is  modeled  as  a  resistor  ( R ). 
The  variations  of  biofluidic  pin  values  due  to  dispersion  effects  are  captured  by  the  model. 

Additionally,  symbols  and  characters  used  in  this  section  are  defined  the  same  as  those  for  the 

mixer,  unless  otherwise  noted.  The  residence  time  of  a  species  band  within  a  separation  channel 
(the  time  for  the  band  to  move  from  the  channel  inlet  to  outlet)  is  given  by 


(21) 


The  calculation  of  changes  in  the  skew  coefficients  and  variance  depends  on  the  specific  ele¬ 
ment  [37]  and  the  inherent  variable  is  the  residence  time  a t  obtained  by  Eq.(21).  For  a  straight  sep¬ 
aration  channel, 


otut-ol  =  Act2  =  2D -A* 

For  a  separation  turn, 
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where  subscrip ts/superscripts  in  and  out  represent  quantities  at  the  inlet  and  outlet  of  the  channel.  In 
Eqs.  (24)  and  (25),  the  “+”  sign  is  assigned  to  the  first  turn  and  any  turn  strengthening  the  skew 
caused  by  the  first;  the  sign  is  assigned  to  any  turn  undoing  the  skew  from  the  first. 

Joule  heating  induced  dispersion  models  were  also  developed  in  this  project  [61],  however, 
are  not  presented  for  the  sake  of  conciseness.  They  are  more  often  than  not  used  to  verily  that  the 
electric  fields  are  not  too  high  to  cause  additional  band  broadening. 

Assuming  a  Gaussian  distribution  of  the  average  concentration  cm  of  the  species  band  at  ele¬ 
ment  terminals  always,  we  can  obtain  the  amplitude  of  the  species  band  by 

KjK  =s\cjI/g:u,  (26) 

For  the  detector  model,  the  variance  change  associated  with  the  detector  path  length  Ldet  is 

given  by  [62] 

Ac2  =  Zjet/l2  (27) 

4,1.3.  Model  Implementation 

To  demonstrate  use  of  the  above  parameterized  models  for  top-down  designs,  we  have 
implemented  the  models  in  the  Verilog-A  analog  hardware  description  language,  and  in  Matlab.  In 
the  Cadence  implementation,  symbol  views  for  each  of  the  elements  are  used  to  compose  a 
schematic  within  Cadence’s  [63]  integrated  circuit  design  framework.  The  Cadence  design 
framework  is  used  to  automatically  netlist  the  complex  topologies  in  the  biofluidic  LoC  schematics, 
and  Spectre  is  used  as  the  simulator.  Similar  tools  from  other  vendors,  or  custom  schematic  entry 
tools  and  solvers  that  can  handle  both  signal  flow  and  Kirchhoffian  networks  could  have  been  also 
used. 

An  important  issue  of  implementing  separation  channel  models  of  turn  geometry  (Eqs.  (24) 
and  (25))  is  the  real-time  detennination  of  the  turn  “sign”.  Providing  this  flexibility  allows  a  single 
turn  symbol  to  be  used  for  constructing  arbitrary  topologies  such  as  a  serpentine,  spiral  or  their 
combination  thereof  as  will  be  shown  later.  To  address  this,  two  sets  of  flags  are  used  in  the  models. 
One  is  the  system  flag  Fs,  stored  as  the  zero -th  component  of  the  skew  coefficients  f  S'[0]  in  Table  1) 
to  record  the  direction  of  the  skew  caused  by  the  first  turn  or  elbow.  The  other  is  the  intrinsic  flag  Fj 
of  individual  elements.  For  example,  Ft  =  1  is  for  turns  or  elbows  involving  clockwise  flow  of 
species  bands;  Ft  =  2  is  for  counter-clockwise  turns  or  elbows.  Since  straight  channels  do  not  incur 
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Table  1. Definition  of  Biofluidic  Pins 


Micromixing 

Bus 

Pins  connected 

Description 

d  [0:29] 

Concentration 

coefficients 

d  [0:9]:  the  1st  sample, 
d  [10:19]:  the  2nd,  d  [20:29]:  the  3rd 

Electrophoretic  Separation 

Bus 

Pins  connected 

Description 

t  [0:2] 

Separation  time 

t[0]  for  the  1st  species,  t[l]  the  2nd,  t[ 2]  the  3rd 

cr  [0:2] 

Variance 

a2-  [0]  for  the  1st  species,  o2  [1]  the  2nd,  cr  [2]  the  3rd 

A  [0:2] 

Amplitude 

A  [0]  for  the  1st  species,  A  [1]  the  2nd,  A  [2]  the  3rd 

S  [0]:  the  direction  of  the  skew  caused  by  the  1st  turn 

S[0:  30] 

Skew  coefficients 

S  [1:10]:  the  1st  species 

5  ri  1:201:  the  2nd.  5  T2 1:301:  the  3rd 

any  skew,  no  flag  is  needed.  During  simulations,  Fs  =  0  (i.e.,  S[0]=0)  is  first  generated  by  the 
injector,  which  is  the  most  upstream  element  of  a  separation  channel  and  hence  initiates  the 
computation  of  the  separation  state.  Then  as  the  species  band  migrates  to  the  first  turn  or  elbow,  Fs 
is  irreversibly  set  to  the  intrinsic  flag  Ft  of  them.  Afterwards,  the  written  Fs  is  compared  with  Ff  of 
each  downstream  element  as  the  band  moves  on.  If  they  are  identical,  a  “+”  sign  is  used  for  the 
element,  otherwise  a  sign.  Fig.  12  shows  the  codes  for  a  turn  involving  clockwise  flow  of 
species  bands  to  implement  this  logic  and  determine  the  sign. 

4.1.4.  Schematic -Based  Simulation 

In  this  section,  we  will  first  describe  the  simulation  procedure,  in  which  the  Kirchhofif’s 
resistor  network  to  predict  electric  current  and  field  as  well  as  the  signal  flow  network  to  evaluate 
biofluidic  state  values  (e.g.,  steady-state  mixing  concentrations  and  transient  electrophoretic  species 
band  shapes)  are  solved  sequentially.  Then,  the  results  of  schematic  simulations  exploring  various 
micromixers  and  separation  microchips  will  be  discussed  and  validated  with  numerical  and  experi¬ 
mental  data. 

Simulation  Description.  Schematic  simulations  for  mixers  and  separation  chips  involve  both  elec¬ 
tric  and  biofluidic  calculations.  For  DC  analysis,  given  the  applied  potential  at  reservoirs,  system 
topology  and  element  dimensions,  nodal  voltages  at  element  tenninals  within  the  entire  system  are 
first  computed  by  Ohm’s  and  Kirchhoff  s  laws  using  the  resistor  models  presented  in  the  last  section. 
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'define  NUM  10 


//  The  number  of  terms  in  the  Fourier  series 


module  Utuml (...); 
parameter  real  sp_D_0  =  500; 
parameter  real  R2  =  1 1 0  ; 
parameter  real  R1  =  100  ; 


//Module  declaration 

//  Diffusivity  of  the  zero-th  species;  unit:  urn  '/s 
//Outer  radius  of  the  U  shape  turn;  unit:  um 
//Inner  radius  of  the  V  shape  turn;  unit:  um 


integer  Fs; 
integer  Fi; 
real  dummyl; 
real  W; 

real  delta_t_0; 
real  pi; 


//  System  flag 
//Intrinsic  flag 
//Intermediate  variable 
//  Channel  width 

//  The  residence  time  of  the  zero-th  species  within  the  turn;  unit:  s 

//  n 


analog  begin 


//Behavioral  description  begins 


Fi  =  1 ;  //  The  present  turn  is  clockwise 

W  =  abs  (R2-R1); 

Fs  =  SKA(sk_in[0]);  // sk  in  is  the  skew  coefficient;  sk_in[0]  stores  the  system  flag 

pi  =  3.1415926; 


if  (Fs  ==  0)  begin 

SKA(sk_out[0])<+(Fi); 
generate  i  (1,  'NUM,  1)  begin 


//system  flag  is  not  set  yet;  no  turns  or  elbows  upstream 
//set  the  system  flag  to  the  intrinsic  flag 
//this  loop  is  to  calculate  skew  coefficient 
dummyl  =  1- exp(-(2*i-1)*(2*i-1)*pi*pi*delta_t_0*sp_D_0/(W*W)); 

Part1_Sk  =  8*pi*W*W*dummy1/(pow((2*i-1),  4)*pow(pi,  4)*delta_t_0*sp_D_0);//fA<;  1”  part  in  Eq.  (13),  use  ‘ 


end 


end  else  if  (Fs  ==  1 )  begin  //Species  flow  in  the  first  turn  is  clockwise 

SKA(sk_out[0])<+SKA(sk_in[0]);  //Convey  the  system  flag 

generate  i  (1,  'NUM,  1)  begin 

dummyl  =  1-  exp(-(2*i-1)*(2*i-1)*pi*pi*delta_t_0*sp_D_0/(W*W)); 

Part1_Sk  =  8*pi*W*W*dummy1/(pow((2*i-1),  4)*pow(pi,  4)*delta_t_0*sp_D_0);  //use  “+” sign 


end  else  begin  //Species  flow  in  the  first  turn  is  counter-clockwise 

SKA(sk_out[0])<+SKA(sk_in[0]);  //  Convey  the  system  flag 

generate  i  (1,  'NUM,  1)  begin 

dummyl  =  1- exp(-(2*i-1)*(2*i-1)*pi*pi*delta_t_0*sp_D_0/(W*W)); 

Part1_Sk  =  -  8*pi*W*W*dummy1/(pow((2*i-1 ),  4)*pow(pi,  4)*delta_t_0*sp_D_0);  //use  “-’’sign 

V 


end 


use  -  sign 


end 

end  // End  of  behavioral  description 

endmodule  // End  of  module 

FIGURE  12.  Verilog-A  description  for  a  180°  turn  involving  clockwise  flow  of  the  species  band.  It 
determines  the  signs  used  in  Eq.  (24),  as  well  as  the  canceling  and  strengthening  effects  of  the  skew. 

The  resulting  nodal  voltages  and  branch  currents  are  in  turn  used  to  calculate  the  electric  field 
strength  ( E )  and  its  direction  within  each  element,  as  well  as  flow  and  splitting  ratios  at  intersections 
(for  mixers).  With  these  results  and  user-provided  sample  properties  ( D  and  //),  the  sample  speed  is 
then  given  by  u  =  juE  .  Next,  values  of  biofluidic  pins  at  the  outlet(s)  of  each  element  (e.g.,  concen¬ 
tration  coefficients  for  micromixers;  arrival  time,  variance,  skew  and  amplitude  for  separation  micro- 
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chips)  are  determined.  The  process  starts  from  the  most  upstream  element,  typically  the  sample 
reservoirs  in  mixers  and  the  injector  in  separation  chips  in  terms  of  the  directional  signal  flow.  As 
such,  both  electric  and  fluidic  information  in  the  entire  system  are  obtained. 

The  mixer  operates  in  steady-state,  while  transient  evolution  is  critical  in  separation 
channels.  Transient  analysis  can  be  also  conducted  for  separation  chips  that  involve  the  species 
band’s  motion  and  broadening.  An  electropherogram  (average  concentration  cm  vs.  time)  can  be 
obtained  at  the  detector,  yielding  an  intuitive  picture  of  separation  resolution  between  species 
bands.  The  transient  analysis  first  calculates  for  the  DC  operating  points  of  the  amplitude  Adet, 

separation  time  tdet  and  variance  <Jdet  of  the  species  band  at  the  detector  as  described  above.  Based 
on  these  points,  the  actual  read-out  time  is  scanned  and  the  average  concentration  output  cm  is 
calculated.  Assuming  the  species  band  does  not  appreciably  spread  out  as  it  passes  through  the 
detector,  cm  is  given  by 


r  _  A  p  2(<&+A°*)  (28) 

Lm  ~  ^dct  ’  ^ 

where  t  is  the  actual  read-out  time  and  A  cr  is  the  variance  growth  associated  with  detection  and 
given  in  Eq.(27). 

Simulation  Results  and  Discussion.  In  this  section,  simulation  examples  of  complex  EK  passive 
mixers  and  electrophoretic  separation  microchips  will  be  presented  to  verily  the  behavioral  models 
for  biofluidic  elements  and  validate  the  modeling  and  simulation  methodology.  Schematic 
simulations  for  micromixers  are  shown  in  Fig.  13  and  Fig.  14  and  Table  2,  and  those  for 
electrophoretic  separation  systems  are  given  in  Fig.  16-Fig.  21. 

Electro  kinetic  micromixers  and  mixing  networks.  Electrokinetic  focusing  [64],  which  first 
appeared  as  an  important  fluid  manipulation  technique  in  EK  LoC  systems,  also  can  be  applied  to 
speed  up  mixing,  especially  in  reaction  kinetics  studies  [65],  Fig.  13  illustrates  an  EK  focusing 
mixer  and  its  system-level  schematic.  In  the  discussion  below,  subscripts  i,  s  and  o  respectively 
denote  the  middle-input,  side  and  output  mixing  channels.  The  cross  intersection  where  sample  a 
(red)  from  the  input  channel  is  focused  by  buffer  or  sample  b  (blue)  from  both  side  channels,  is 
modeled  as  two  concatenated  converging  intersections.  The  flow  ratio  (the  ratio  of  the  flow  rate  of 

the  middle-input  stream  to  the  total  flow  rate)  of  sample  a  is  s  =  /  /(2/  +  /,.). 
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a 


b 


FIGURE  13.  (a)  An  electrokinetie  focusing  micromixer.  Sample  a,  flowing  from  the  top  input  channel  to  the 
intersection,  is  pinched  by  sample  b  (or  buffer)  from  both  side  channels.  Then  samples  mix  in  the  bottom 
mixing  channel,  (b)  Its  hierarchical  schematic  representation.  The  triple  input  and  one  output  cross 
intersection  is  modeled  as  a  cascade  connection  of  two  converging  intersections. 


Fig.  14(a)  shows  numerical  and  schematic  simulation  results  of  sample  a  concentration  pro¬ 


files  at  the  mixing  channel  outlet  for  two  flow  ratios  s  =  0.1  and  s  =  1/3 .  During  simulations,  reser¬ 
voir  voltages  (</>j  and  (f>s)  are  chosen  to  vary  s  while  keeping  E  (143  V/cm)  and  the  sample  residence 


time  fixed  in  the  mixing  channel.  Excellent  agreement  between  numerical  and  schematic  simulation 


results  is  found  with  the  worst-case  relative  error  of  3%  at  s  =  0.1.  The  results  are  also  compared 


with  those  from  a  T-type  mixer  that  has  the  same  electrical  field  (in  the  mixing  channel),  channel 


0.0  0.2  0.4  0.6  0.8  1.0 


Norm,  widthwise  position  T) 

(a) 


Axial  position  z  (cm) 
(b) 


FIGURE  14.  (a)  Schematic  simulation  results  (lines)  compared  with  numerical  data  (symbols)  on 
concentration  profiles  c  (sample  a)  along  the  normalized  channel  width  77  for  the  electrokinetie  focusing  and 
T-type  mixers.  In  contrast  to  the  T-type  mixer,  the  focusing  mixer  considerably  improves  sample  homogeneity 
due  to  the  reduced  diffusion  distance  between  samples;  (b)  Schematic  simulation  results  on  variation  of 
mixing  residual  Q  along  axial  channel  length  (data  points  are  connected  by  lines  to  guide  the  eye)  for  the 
electrokinetie  focusing  mixer  involving  different  stream  width  s.  A  smaller  stream  width  (e.g.,  s  =0.1 )  yields 
a  lower  initial  mixing  residual  and  a  more  uniform  concentration  profile  at  the  end  of  the  mixing  channel. 
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Table  2.  Comparison  of  schematic  simulation  results  (sche)  with  numerical  (num)  and  experimental  (exp)  data  on 
sample  concentrations  in  analysis  channels  of  serial  and  parallel  mixing  networks  [60] 


Serial  Mixing  Network 

Parallel  Mixing  Network 

Complete  Mixing 

Partial  Mixing 

Complete  Mixing 

channel 

c  (sche) 

c  (exp) 

c  (num) 

c  (sche) 

c  (num) 

channel 

c  (sche) 

c(exp)  c 

•  (num) 

Ai 

1 

1 

1 

1 

1 

Al 

0 

0 

0 

A2 

0.37 

0.36 

0.378 

0.48 

0.496 

A2 

0.83 

0.84 

0.832 

A3 

0.22 

0.21 

0.224 

0.187 

0.187 

A3 

0.68 

0.67 

0.674 

A4 

0.125 

0.13 

0.133 

0.081 

0.0815 

A4 

0.52 

0.51 

0.523 

A5 

0.052 

0.059 

0.0628 

0.029 

0.0315 

A5 

0.35 

0.36 

0.354 

A6 

0.17 

0.19 

0.168 

A7 

1 

1 

1 

length  and  width  as  the  focusing  mixer.  The  focusing  mixer  considerably  improves  sample  homoge¬ 
neity,  which  can  be  attributed  to  the  reduced  diffusion  distance  between  samples.  That  is,  the  axial 
centerline  of  the  mixing  channel  in  the  focusing-mixer  can  be  treated  as  an  impermeable  wall  due  to 
symmetry,  hence  the  diffusion  distance  is  only  one -half  of  that  of  the  T-type  mixer.  Also,  a  smaller 
stream  width  (e.g.,  s  =  0.1)  yields  a  more  uniform  concentration  profile  at  the  end  of  the  mixing 
channel. 

To  gain  the  insight  of  the  influence  of  the  stream  width  on  mixing  performance,  an  index  of 
mixing  residual,  Q  =  \'\c(ri)-cavg\dri ,  is  introduced  in  Fig.  14(b)  to  characterize  the  non-unifonnity 
of  concentration  profiles,  where  c(rj)  and  cavg  are  the  normalized  concentration  profile  and  width- 
averaged  concentration  respectively  at  the  detection  spot.  At  the  channel  inlet  (z  =  o ),  mixing 
residual  Q  strongly  depends  on  5.  Asymmetric  incoming  streams  yield  a  lower  Q  value  (e.g., 
Q  =  0.18  at  =  o.l  compared  with  0  =  0.44  at  .?  =  l/3)  and  a  more  unifonn  initial  profile.  Along  the 
channel,  Q  initially  drops  rapidly  and  then  becomes  saturated  because  the  improved  sample  mixing 
reduces  the  concentration  gradient  and  the  driving  force  for  further  mixing.  Thus,  a  tradeoff 
between  Q  and  mixer  size  can  be  captured  by  the  behavioral  models  presented  in  this  paper  to 
achieve  designs  of  both  effective  and  efficient  micromixers. 

These  parameterized  and  reusable  behavioral  models  are  well  suited  to  study  complex  mix¬ 
ing  networks  [60],  in  which  an  array  of  sample  concentrations  can  be  obtained  at  multiple  analysis 
channels  by  geometrically  duplicating  units  with  a  single  constant  voltage  applied  at  all  reservoirs. 

Table  2  shows  the  comparison  of  schematic  simulation  results  with  experimental  and 
numerical  data  on  sample  (rhodamine  B)  concentrations  in  analysis  channels  A  |-A5  in  the  serial 
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FIGURE  15.  An  electrokinetic  serial  mixing  network  [60]  and  its  hierarchical  schematic  representation.  The 
network  consists  of  reservoirs,  mixing  channels,  T-  and  cross-intersections.  Sample  and  buffer  is  released 
and  collected  by  the  reservoirs.  In  the  composable  approach,  the  serial  mixing  network  is  represented  as  a 
collection  of  interconnected  mixing  elements  composed  of  microchannels,  converging  intersections  and 
diverging  intersections. 

mixing  network  (Fig.  15).  Both  complete  and  partial  mixing  cases  are  investigated.  When  a  voltage 
of  0.4  kV  is  applied  at  the  sample  and  buffer  reservoirs  with  the  waste  reservoirs  grounded,  sample 
mixing  in  channels  S2-S5  is  width-wisely  complete.  Excellent  agreement  among  the  schematic 
simulation,  numerical  analysis  and  experimental  results  is  found  (with  an  average  error  smaller  than 
6%  relative  to  experiments).  In  contrast  to  the  electric  resistor-based  models  [22]  [32]  [60]  that  take 
advantage  of  the  analogy  between  EK  flow  and  electric  current  and  require  post-calculations  of 
concentrations  from  current  distributions  in  the  network,  our  behavioral  models  directly  yield  the 
concentration  value  in  each  analysis  channel.  In  addition  to  complete  mixing,  partial  mixing  case  is 
also  schematically  simulated.  A  voltage  of  1.6  kV,  as  used  in  the  experiments  in  [60],  is  applied  at 
the  sample  and  buffer  reservoirs  with  the  waste  grounded,  which  increases  the  EK  velocity  and  then 
decreases  the  residence  time  of  the  sample  by  four  folds  in  channels  S2-S5.  Thereby,  the  mixing  in 
channels  S2-S5  is  width-wisely  incomplete,  and  the  amount  of  sample  shunted  to  channels  A|-A5 
depends  on  not  only  the  electric  currents  in  the  network  but  also  the  sample  concentration  profiles  at 
the  exits  of  channels  S2-S5,  which  violates  the  assumption  for  the  analogy  between  EK  flow  and 
electric  currents  and  hence  the  resistor-based  modeling  becomes  invalid.  However,  it  can  be  readily 
simulated  by  our  behavioral  models.  In  the  schematic,  the  cross-intersection  is  modeled  as  a 
combination  of  the  converging  and  diverging  intersections,  in  which  the  sample  concentration 
profiles  of  the  incoming  and  outgoing  streams  are  accurately  captured.  Results  from  schematic 
simulations  are  compared  with  numerical  data  in  Table  2  (a  comparison  to  experimental  data  is  not 
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allowed  due  to  a  lack  of  knowledge  of  sample  properties.  Hence  a  difiusivity  of  £>=3x10  nr/s 

and  an  EK  mobility  of  //  =2.0x10  m  /Vs  are  assumed  in  numerical  simulations).  Very  good 
agreement  can  be  observed  with  an  average  error  of  4%.  At  the  cross-intersection  following  channel 
S2,  the  amount  of  sample  shunted  to  A2  is  more  than  that  predicted  by  the  complete -mixing  case 
due  to  the  non-uniform  sample  profiles  at  the  intersection’s  inlet.  Consequently,  concentrations  in 
channels  A3-A5  show  the  lower  values,  which  agrees  with  experimental  observations  [60]. 
Netlisting  and  schematic  simulation  of  this  example  take  20  seconds  on  a  multi-user,  2-CPU  1-GHz 
Sun  Fire  280  processors  with  4  GB  RAM  for  the  first  time  simulation,  and  less  than  a  second  for 
subsequent  iterations,  leading  to  a  1 000-20, OOOx  speedup. 

In  additional  to  the  serial  mixing  network,  the  parallel  mixing  network  [60]  can  be  hierarchi¬ 
cally  represented  and  simulated  in  a  similar  fashion  and  excellent  agreement  among  schematic  sim¬ 
ulations  results,  numerical  analysis  and  experimental  data  (with  an  average  error  of  3.6%  relative  to 
experiments)  is  also  found  [59]. 

Electrophretic  separation  microchips.  Schematic  simulation  results  for  electrophoretic  separation 
microchips  are  shown  in  Fig.  16-Fig.  21.  In  Fig.  16  and  Fig.  17,  a  serpentine  electrophoresis 
column  of  two  complementary  turns  is  used  to  separate  an  analyte  comprised  of  two  species  a 
(Z)=3.12xl O’10  nr/s,  //=1.2xl0"8  irr/sV)  and  b  (£>=2.72xl0_1°  nr/s,  //=l.lxl0"8  nr/sV)  with 
E  =  600  V/cm.  Experimental  data  [66]  on  variance  vs.  time  of  species  a  are  compared  with  DC 
schematic  simulations  in  Fig.  16,  showing  excellent  agreement  with  the  worst-case  relative  error  of 
only  5%.  Again,  netlisting  and  DC  simulation  for  this  example  take  20  seconds  for  the  first  time  and 
less  than  a  second  for  subsequent  iterations,  leading  to  a  500- 10, OOOx  speedup  (higher  speedup  can 
be  obtained  for  a  more  complex  chip  topology  or  a  less  diffusive  species  as  shown  in  Fig.  20).  The 
first  turn  skews  the  species  band  and  accordingly  incurs  abrupt  increase  in  variance.  During  the 
species  band’s  migration  in  the  long  inter-turn  straight  channel,  the  transverse  diffusion  smears  out 
most  of  the  skew  and  presents  a  nearly  unifonn  band  before  the  second  turn.  The  second  turn  then 
distorts  the  band  again  in  the  opposite  direction,  leading  to  another  turn-induced  variance  that  is 
equal  to  the  one  caused  by  the  first  turn.  Fig.  17  shows  electropherograms  of  both  species  from 
three  detectors.  The  spacing  between  concentration  peaks  of  species  a  and  b  increases  as  they 
migrate  through  channels,  but  due  to  the  band-broadening  effect,  the  amplitude  decreases 
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FIGURE  16.  Comparison  of  experimental  data  [66]  with  DC  schematic  simulation  on  variance  a2  vs. 
separation  time  t  of  species  a  in  a  serpentine  electrophoretic  separation  microchip  of  two  complementary 
turns.  The  grey  bars  represent  the  residence  time  of  the  sample  within  the  turns.  The  first  turn  skews  the 
species  band  and  accordingly  incurs  abrupt  increase  in  variance.  The  transverse  diffusion  in  the  inter-turn 
straight  channel  smears  out  most  of  the  skew  and  presents  a  nearly  uniform  band  before  the  second  turn  (see 
the  inset  of  numerical  simulation  plot).  The  second  turn  then  distorts  the  band  again  in  the  opposite 
direction,  leading  to  another  turn-induced  variance. 


FIGURE  17.  Transient  analysis  simulates  the  electropherograms  output  from  three  detectors,  which  are 
respectively  arranged  before  the  first  turn  (top  trace),  in  the  middle  of  the  inter-turn  straight  channel 
(middle  trace)  and  after  the  second  turn  (bottom  trace).  Attributed  to  the  difference  in  electrokinetic 
mobility,  the  spacing  between  concentration  peaks  of  species  a  and  b  increases  as  they  migrate  through 
channels.  The  dispersion  effect  leads  to  the  continuous  decreases  in  the  band  amplitude. 
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FIGURE  18.  (a)  A  spiral  electrophoretic  separation  microchip  |20].  It  consists  of  five  turns  with 
continuously  decreased  radius  (1.9, 1.8, 1.7, 1.6  and  0.8  cm).  Within  them,  species  Dichlorofluoroscein  flows 
in  the  same  direction  (clockwise),  (b)  Its  hierarchical  schematic  representation. 


consecutively. 

In  Fig.  18,  the  dispersion  of  Dichlorofluoroscein  in  a  complex  spiral  separation  microchip  of 
five  turns  is  simulated  and  compared  with  experimental  results  [20].  Spiral  channels  differ  from  the 
serpentine  in  that  the  skew  and  variance  always  increase  with  the  turn  number,  as  the  band  skew  in 
all  turns  has  the  same  sense  and  does  not  cancel.  A  scalar  index  of  plate  number  Ns  to  characterize 

the  resolving  power  of  the  electrophoresis  chip  is  defined  Ns  =l}toJa2 ,  where  Ltot  is  the  total 
separation  length  from  injector  to  the  detector.  The  higher  the  plate  number,  the  better  separation 
performance  achieved  by  the  chip.  The  linear  growth  of  the  plate  number  with  electric  field  implies 
that  molecular  diffusion  is  the  major  dispersion  source  in  such  a  system  Fig.  19),  as  molecular 
diffusion  decreases  linearly  as  electric  field  increases.  The  worst-case  relative  error  of  12%  is 
considered  acceptably  small  considering  the  uncertainties  in  the  knowledge  of  species  diffusivity 
[20]. 

Fig.  20  illustrates  a  hybrid  electrophoretic  separation  microchip  [67]  and  its  schematic 
representation  including  both  spiral  and  serpentine  channels.  Due  to  the  difficulty  of  accounting  for 
the  coexisting  skew  canceling  and  strengthening  effects  in  such  a  topology,  it  has  not  been 
effectively  investigated  since  it  was  proposed  [67].  Fig.  21  shows  schematic  simulation  result  on  the 
variance  of  a  species  band  vs.  time  in  such  a  chip,  as  well  as  its  comparison  with  numerical  data.  A 

low  species  diffusivity  0f  £>-1x10  m  /s  is  chosen  to  analyze  the  highly  convective  dispersion 
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FIGURE  19.  Comparison  of  schematic  simulation  results  to  the  experimental  data  on  theoretical  plate 
number  Ns  vs.  electric  field  E.  Right  axis  shows  the  relative  error  between  simulation  and  experiments.  The 
linear  growth  of  the  plate  number  Ns  with  electric  field  E  implies  that  molecular  diffusion  is  the  major 
dispersion  source  in  such  a  system. 


that  has  not  been  considered  by  the  previous  example  in  Fig.  16  (other  properties  and  parameters 
are  the  same  as  those  of  sample  a  in  Fig.  16).  Highly  convective  dispersion  is  practically  important 
for  microchip  electrophoresis  of  the  species  with  low  diffusivity,  such  as  the  separation  of  DNA  in  a 


FIGURE  20.  (a)  A  hybrid  electrophoretic  separation  microchip.  It  consists  of  both  spiral  and  serpentine 
channels.  Species  flows  in  the  clockwise  direction  in  both  turns  Tj  and  T2  (spiral  topology),  thereby  T2 
strengthens  the  sharp  skew  generated  by  Tj.  The  skew  almost  persists  through  the  inter-turn  straight 
channel  between  T2  and  T3  and  is  significantly  cancelled  out  by  T3  (serpentine  topology),  (b)  Its 
hierarchical  schematic  representation. 
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Separation  time  t  (s) 

FIGURE  21.  Comparison  of  numerical  data  with  DC  schematic  simulation  on  variance  vs.  separation  time 
in  the  hybrid  electrophoretic  separation  microchip.  Very  sharp  skew  (see  Fig.  20)  is  generated  and  the 
variance  accumulates  after  turns  Tj  and  T2  due  to  their  spiral  topology.  The  skew  almost  persists  through 
the  inter-turn  straight  channel  between  T2  and  T3  and  is  significantly  cancelled  out  by  T3  attributed  to  their 
serpentine  topology,  which  as  a  result  yields  a  drastic  variance  drop  after  T3. 

gel  or  sieving  matrix  [66] [68].  It  is  shown  in  Fig.  20  that  since  species  flows  in  the  clockwise 
direction  in  both  turns  T  |  and  To  (spiral  topology),  To  strengthens  the  sharp  skew  generated  by  T1? 
leading  to  a  more  skewed  band  and  a  higher  variance.  Due  to  the  small  species  diffusivity,  the  skew 
almost  persists  through  the  inter- turn  straight  channel  between  To  and  T3  and  is  significantly 
cancelled  out  by  T3,  which  as  a  result  yields  a  drastic  variance  drop  in  T3  (serpentine  topology). 
However,  the  skewed  band  after  T3  is  overly  corrected  by  T4  and  a  counter-skew  is  shown 
afterward.  Excellent  agreement  between  the  schematic  and  numerical  simulation  results  with  1% 
relative  error  and  tremendous  speedup  up  to  400,000x  have  been  achieved  in  Fig.  21.  This  is  the 
first  time  that  the  highly  convective  dispersion  in  the  hybrid  electrophoresis  microchip  at  this 
complexity  level  has  been  accurately  and  efficiently  simulated  by  analytical  models. 

4.1.5.  Injector  Models 

Injector  Basic  Operation  and  Topologies.  The  injector  prepares  the  sample  from  the  synthesis 
stage  of  the  LoC  for  measurement  in  the  separation  stage.  This  is  done  very  simply  by  the  intersec¬ 
tion  of  one  or  more  channels  as  shown  in  Fig.  22.  It  should  be  noted  that  injectors  are  classified  not 
only  by  topology,  but  also  by  the  field  arrangements  used  to  shape  the  fluids  inside— for  example  the 
cross  and  gated-cross.  The  earliest  on-chip  injectors  were  the  simple  tee  [69][70].  The  tee  operates 
in  two  stages  by  first  flowing  analyte  past  the  empty  separation  channel.  Once  enough  material  has 
passed  the  separation  channel,  the  flow  is  reversed  by  modifying  the  potential  on  the  separation 
channel  to  create  an  electric  field,  capturing  a  plug  of  material.  This  method  has  poor  control  over 
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the  size  of  the  material  plug,  [71],  and  has  mostly  be  replaced  by  the  cross,  [72],  and  double-tee, 
[73].  The  cross  and  double-tee  also  operate  in  two  stages  and  they  operate  on  similar  principles;  the 
double-tee  simply  provides  a  larger  injected  band  due  to  the  channel  offset.  In  both  injectors,  mate¬ 
rial  is  first  moved  into  the  channel  intersection,  called  the  injection  chamber.  While  material  is  mov¬ 
ing  into  the  chamber,  electric  fields  from  the  side  channels  (accessory  fields)  constrict  the  flow  to 
provide  a  thinner  band  to  the  separation  channel.  In  the  second  stage,  the  fields  change  direction  to 
flush  the  thin  column  of  material  into  the  separation  channel.  Simultaneously,  the  side  channels 
(now  the  top  and  bottom  channels)  have  fields  directed  away  from  the  injection  chamber  to  prevent 
additional  material  from  flooding  into  the  separation  channel,  guaranteeing  a  small  clean  plug 
[62] [74] [75].  The  gated-cross  has  the  same  topology  as  the  regular  cross,  but  operates  with  signifi¬ 
cantly  different  electric  fields,  [62][76][77][78][79].  The  gated-cross  is  desirable  when  repeated 
injections  must  be  performed  in  rapid  succession.  The  injector  operates  by  running  the  material  at  a 
90  degree  angle  from  the  input  channel  using  an  opposing  field  to  fonn  a  gate.  The  gate  is  released 
for  a  given  time  to  obtain  a  plug  of  the  desired  size  and  then  closed.  When  the  gate  is  closed  the 
opposing  field  flushes  the  plug  into  the  separation  channel.  Meanwhile  the  original  material  contin¬ 
ues  to  flow  on  the  other  side  of  the  gate  and  is  immediately  ready  for  another  injection.  The  double- 
cross  injector  is  similar  to  the  regular  cross,  but  includes  an  extra  set  of  accessory  channels  to  create 
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FIGURE  22.  Basic  operating  stages  for  various  injector  topologies.  The  most  common  topologies  are  the 
cross,  double-tee,  and  gated  cross.  The  earliest  on-chip  injectors  were  the  tee. 


38 


an  extra  thin  and  uniform  stream  for  plug  extraction,  [80].  The  additional  complexity  of  controlling 
six  channels  simultaneously  is  typically  not  worth  the  incremental  improvement  over  the  regular 
cross,  and  thus  this  topology  has  not  seen  wide  acceptance.  Other  lesser-used  topologies  and  control 
methods  have  been  explored  for  specific  applications  such  as  the  multi-legged  fork  [81],  which  pro¬ 
vides  discrete  volumes  based  on  the  number  of  legs  fabricated  and  the  distance  of  the  inflow  and 
outflow  legs,  and  the  gated-double-tee  and  triple-tee  [82],  which  provide  modifications  to  the  regu¬ 
lar  double-tee. 

The  goal  of  all  of  these  injection  methods  is  to  increase  the  efficiency  of  the  separation  chan¬ 
nel.  The  injector  should  provide  a  small  well-defined  plug  of  material  with  maximum  concentration 
for  ease  of  separation.  The  smaller  the  plug  is,  the  shorter  the  distance  each  species’  peak  center 
must  travel  to  be  resolvable,  thus  shorter  channels  can  be  used,  allowing  room  for  more  complex 
system  integrations.  An  example  of  this  is  efficiency  shown  in  Fig.  23.  The  parameters  of  the  simu¬ 
lation  are  shown  in  the  caption  of  the  figure,  and  the  results  show  that  a  wider  band  requires  a 
longer  separation  channel  in  order  to  get  the  same  resolution. 


(a) 


(b) 


I  I 


FIGURE  23.  Two  species  are  combined  and  a  band  of  material  representative  of  an  injector  formed  plug  is 
placed  at  the  beginning  of  the  channel.  Species  one  and  two  have  a  diffusivity  of  lxlO'10  m2/s.  Species  one 
has  a  mobility  of  5.25xl0'7  m2/Vs,  whereas  species  two  has  a  mobility  of  5.5x10-7  m2/Vs.  After  traveling  in 
a  50000  V/m  electric  field  for  1.2mm,  the  separation  of  the  two  species  bands  is  measured.  In  the  first 
example  (a),  the  plug  has  a  width  of  lOum.  In  the  second  example  (b),  the  plug  has  a  width  of  50um.  Even 
though  the  plug  in  (b)  contains  more  material,  it  is  not  as  easy  to  resolve  as  the  smaller  plug  in  (a). 
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Injector  Modeling  Parameters.  The  parameters  can  be  categorized  into  one  of  three  categories: 
(1)  Those  describing  the  fluidic  properties  of  the  analyte,  (2)  those  describing  the  applied  electric 
fields,  and  (3)  those  describing  the  geometry  of  the  given  injector's  topology.  The  fluidic  properties 
include  the  physicochemical  attributes  of  the  buffer  fluid  and  the  transported  chemical  species  such 
as  electrokinetic  mobilities  and  diffusivities.  The  applied  electric  fields  provide  the  stimulus  and 
driving  forces  for  the  fluids.  These  fields  are  characterized  by  their  direction  and  magnitude  in  each 
of  the  channels.  In  general,  the  conservation  of  current  and  ohmic  conduction  are  inherent  to  the 
problem  so  although  there  may  be  N  channels,  there  are  only  N- 1  independent  field  quantities.  For 
each  stage  of  operation  the  field  directions  or  magnitudes  may  vary,  so  for  M  stages  of  operation 
there  are  M(N- 1)  independent  field  parameters.  The  relevant  geometric  parameters  for  injectors  are 
the  channel  widths  and  lengths,  where  transport  in  the  depth  direction  is  typically  ignored.  The 
widths  of  the  injector  channels  near  their  intersection  point  are  almost  always  the  same,  so  knowing 
one  width  provides  information  about  all  channel  widths.  The  effects  of  unequal  channel  widths  at 
the  intersection  has  been  explored  [83].  Unequal  channel  widths  has  the  ability  to  enable  injections 
with  fewer  controlling  electrodes,  but  is  not  commonly  implemented  due  in  part  to  the  difficulty  of 
filling  such  channels  with  buffer  when  dry.  The  unequal  channel  resistances  make  the  pressure 
driven  filling  more  difficult. 

Finally,  unless  noted  otherwise  the  injector  is  treated  as  a  steady  state  device.  Each  stage  of 
operation  is  run  until  stable  equilibrium  is  reached.  For  the  cross  and  double -tee  loading  stages,  this 
equilibrium  is  achieved  in  all  cases,  except  when  the  side  channels  are  left  floating  (no  accessory 
electric  fields),  permitting  diffusion  into  the  side  channels  that  floods  the  chip.  In  practice,  the  injec¬ 
tor  is  only  loaded  for  a  limited  time  to  prevent  flooding  the  side  channel  with  diffusion  in  the 
absence  of  accessory  fields.  In  general,  this  time  is  long  enough  to  reach  a  stable  steady  state  for  the 
instances  when  there  are  accessory  electric  fields  in  the  side  channels.  In  this  work  the  following 
time  is  used: 


t  =  8.4^—  + 

L  W, 


5Z  +  ^3 

ul2av  u, 


(29) 


where  Z;  is  the  length  traveled  from  the  loading  stage  input  and  waste  port  to  the  channel  intersec¬ 
tion,  Uj  is  the  velocity  in  the  input  channel,  A  is  the  offset  length  of  the  channels— zero  for  the  cross, 
U iiav  is  the  average  velocity  of  the  input  channel  and  the  waste  channel,  Z?  is  the  length  from  the 
channel  intersection  to  the  waste  port,  and  U3  is  the  velocity  in  the  waste  channel.  The  factors  of  8.4 
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and  5  are  empirically  chosen  to  ensure  sufficient  time  for  steady  state  to  be  achieve  when  accessory 
fields  are  implemented. 

The  loading  stage  is  characterized  by  a  time  that  allows  the  band  to  completely  exit  the 
injection  chamber: 


_  L d 

d  -  Tj, 


(30) 


where  Lj  is  a  characteristic  length  for  the  band  to  travel  to  exit  the  injection  chamber  and  Uci  is  the 
velocity  of  the  band  in  the  separation  channel.  Beyond  this  point,  the  effects  of  the  separation  channel 
determine  the  dispersion  of  the  band. 


Injector  Leakage.  For  proper  operation  or  simulation  all  of  the  parameters  and  their  appropriate 
ranges  must  also  be  defined.  This  defines  the  design  space  for  the  injector  subsystem.  An  injector  is 
only  properly  operating  for  certain  values  of  the  electric  field,  and  the  configuration  of  electric 
fields  contributes  to  the  classification  of  an  injector.  For  example,  if  the  input  electric  field  reverses 
direction  and  prevents  the  injection  chamber  from  filling  with  chemical  species  then  the  device  will 
not  have  the  material  necessary  to  form  a  plug  for  separation.  Another  important  implication  is  that 
the  fields  must  be  applied  in  the  proper  ratio  for  the  dispensing  stage  to  prevent  flooding  the  separa¬ 
tion  channel  with  material  from  the  upstream  mixer  or  reservoir.  This  is  achieved  by  requiring  the 
field  in  side  channels  during  dispensing  to  be  larger  than  some  multiple,  a,  of  the  driving  field  in  the 
separation  channel.  This  is  depicted  in  the  cartoon  of  Fig.  24,  and  the  factor,  a,  depends  on  the 
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FIGURE  24.  Proper  operation  of  an  injector  requires  the  identification  of  control  parameters  and  their 
applicable  ranges.  Here  identifying  the  fields  (arrows)  is  the  first  step,  but  then  the  magnitude  and  ratios 
must  be  additionally  constrained  for  proper  operation  ( E1L  >  0,  E1D/E2D  >  a  ,  where  E1D=E3D) 
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Peclet  number  for  this  dispensing  stage,  since  more  diffusion  will  be  conducive  to  more  leakage. 
Applying  the  Buckingham-Pi  theorem  identifies  the  essential  non-dimensional  parameters  to  be  a 
field  ratio  in  the  loading  stage  (L)  and  dispensing  stage  (D),  sL  =  E2L/E3L,  eD  =  E1D/E2D,  and  the 
Peclet  numbers  in  the  loading  and  dispensing  stage,  PeL  =  Ujw/tc,  PeD  =  U2w/k.  Where  k  is  the  dif- 
fusivity,  Uj  is  the  velocity  in  channel  i,  and  w  is  the  channel  width. 

To  measure  a,  simulations  were  done  for  a  cross  injector.  To  quantify  whether  there  was  sig¬ 
nificant  leakage,  the  concentration  was  measured  at  the  immediate  edge  of  the  injection  chamber 
adjacent  to  the  separation  channel,  after  the  band  had  traveled  six  channel  widths  downstream.  If 
the  ratio  of  the  peak  of  the  injected  band  to  the  value  at  the  edge  of  the  injection  chamber,  exceeds 
the  ratio  of  the  peak  of  a  Gaussian  distribution  to  its  value  at  3a,  leakage  is  defined  as  excessive. 
The  following  relation  defines  the  region  of  acceptable  leakage  operation: 


r 

max  y 

C  “ 

min 


(31) 


where  Cmax  is  the  peak  of  the  injected  band  and  Cmin  is  the  concentration  at  the  exit  of  the  injection 
chamber. 


The  leakage  obviously  depends  on  the  ratio  of  the  accessory  fields  in  the  side  channels  to  the 
driving  field  in  the  separation  channel  and  the  Peclet  number  for  the  dispensing  stage.  To  test  for 
dependence  on  the  loading  stage  parameters,  the  field  ratios  and  Peclet  number  for  the  loading  stage 
were  varied  as  well.  The  results  are  displayed  in  Fig.  25  and  show  a  slight  dependence  on  the  load¬ 
ing  stage  field  ratios.  The  more  constricted  the  band  is  during  the  loading  stage,  the  less  leakage  is 
likely  in  the  dispensing  stage.  The  worst  case  (low  loading  field  ratios  and  low  loading  Peclet  num¬ 
ber)  is  used  to  create  a  bound  to  restrict  the  design  space  of  the  cross  and  double-tee  injectors: 


1  ( PeD\ 

S°  >  ^IfT65  \849o)  (32) 

Injector  Output  Plug  Gaussian  Representation.  The  size  and  shape  of  the  plug  of  material 
injected  into  the  separation  channel  is  the  modeled  output  of  the  injector.  The  numerical  results  sim¬ 
ulated  at  all  of  the  points  defined  by  the  injector  design  space  contain  all  of  the  infonnation  regard¬ 
ing  the  plug.  To  efficiently  represent  the  plug,  the  results  of  the  numerical  simulation  are  processed 
to  extract  the  equivalent  Gaussian  peak-height  and  variance.  The  injector  model  has  the  following 
form: 
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where  n  are  the  input  parameters  for  the  specific  injector  device  being  modeled,  cr  is  the  variance 
of  the  widthwise  average  concentration  from  the  numerically  simulated  injector,  y\r  is  the  channel 
width  normalizing  factor,  C0  is  the  normalizing  concentration  —  typically  the  output  concentration 
of  the  upstream  mixer,  and  Cp  is  the  peak-height  of  an  equivalent  Gaussian  for  the  transversely  aver¬ 
aged  concentration  of  the  numerically  simulated  injector  such  that  the  effective  Gaussian  and  the 
actual  plug  have  the  same  mass, 
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m 
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where  m  is  the  mass  of  the  band,  and  cr  is  the  band  variance. 


(34) 


This  is  an  efficient  representation  of  the  band  since  the  effects  of  diffusion  will  always  act 
asymptotically  to  make  any  band  Gaussian.  However,  since  diffusion  takes  time  to  homogenize  the 
band,  there  is  a  limit  on  the  applicability  of  this  method.  This  bound  can  be  identified  based  on  sev¬ 
eral  time  constants  associated  with  the  downstream  separation  channel.  The  effective  Gaussian  rep¬ 
resentation  will  be  an  accurate  representation  as  long  as  the  detector  is  far  enough  downstream  to 


H-O.Pe,  =0  =  0,  Pe,  =  1000 


%  % 


FIGURE  25.  Leakage  results  for  the  cross  and  double-tee  injectors.  The  values  of  the  loading  stage 
parameters  are  fixed,  but  different  in  each  of  the  four  plots.  The  blue  region  defines  infeasible  operation, 
and  the  red  region  defines  feasible  operation  for  the  dispensing  stage  field  ratios,  sD  and  Peclet  number 
PeD- 
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L  =  0  :  r=  °o  L  =  4.7  mm  :  r=  2 


FIGURE  26.  Comparison  of  actual  injection  plug  to  its  effective  Gaussian  representation.  The  actual 
output  of  a  double-tee  injector  is  in  the  top  channels,  the  effective  Gaussian  model  is  shown  in  the  bottom 
channels.  The  band  on  the  right  has  traveled  4.7mm,  which  is  when  r=  2.  For  these  simulation  n1D  =  186, 
n2D  =  186,  Ji3D  =  1/8,  and  7i4D  =  0.57,  ji5D  =  2  (definitions  can  be  found  in  Table  3. 

allow  diffusion  to  work  and  before  there  are  any  additional  sources  of  dispersion  such  as  turns  in  the 

channel.  In  other  words,  the  time  it  takes  particles  to  diffuse  across  the  width  of  the  channel,  w, 


t 


K 


must  be  shorter  than  the  time  it  takes  particles  to  reach  the  first  channel  bend  or  detector,  L, 


(35) 
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Using  these  times  and  requiring  tc>tv  the  following  bound  can  be  used  to  ensure  that  the 


band  is  accurately  represented  by  a  Gaussian  before  being  detected  at  the  end  of  the  separation 
channel: 
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(37) 


where  Pew  is  the  Peclet  number  based  on  the  channel  width.  One  implication  of  this  bound  is  that  it 
not  accurate  to  query  the  detailed  shape  of  the  plug  immediately  after  the  injector,  since  the  actual 
plug  has  not  had  a  chance  to  homogenize  into  its  ultimate  Gaussian  shape. 


As  an  example  of  the  effectiveness  of  this  representation,  Fig.  26  shows  the  result  of  a  band 
created  from  a  double -tee  injector  as  it  travels  down  a  straight  separation  channel.  Immediately 
after  the  injector  the  differences  in  the  plug  are  apparent,  but  as  the  actual  band  travels  through  the 
channel  it  becomes  Gaussian  due  to  the  effects  of  diffusion  and  is  then  indistinguishable  from  the 
effective  Gaussian  plug. 


Cross  Topology.  The  standard  operation  for  a  cross  injector  is  shown  in  Fig.  27.  It  is  a  two  stage 
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Table  3. Cross  Injector  Simulation  Parameters 


Simulation  Parameters 

System  Non-Dimensional 

Non-Dimensional  Parame- 

Parameters 

ter  Constraints 

Incut 

Input 

Input 

H  -  Electrokinetic  mobility 

nlc  =  E2L,/ElL 

n3c,  n4c  e  [10,  5000] 

K  -Diffusion  coefficient 

k2c=  EmJEzn 

nlc,  n2c  e  [0.01,  5] 

Eil,  E2l,  E3l,  E4L-  Load¬ 
ing  electric  field 

"3c  =  JuE21w/K 

3 

cq 

II 

si 

E id>  ^2D’  E2d,  E4D  -  Dis¬ 
pensing  electric  field 

ti4c=  /jE2dw/K 

E1D-E3D 

w  -  Channel  width 

n2c>  1/16.65  log(8490/7t4c) 

C0  -  Analyte  concentration 

Output 

Output 

Cp  -  Equivalent  Gaussian 
Peak 

Cp=Cp/C0 

a2  -  Equivalent  Gaussian 
Variance 

a2  =  a2/w2 

device  where  the  material  from  the  upstream  mixer  is  first  loaded  into  the  intersecting  cross  channel 
region.  This  is  done  using  one  configuration  for  the  electric  fields.  Once  this  region  has  been  filled, 
the  electric  fields  are  modified  to  drive  a  plug  of  fluid  into  the  separation  channel.  Of  the  fluidic, 
electric,  and  geometric  parameters  mentioned  above,  the  ones  relevant  to  the  cross  injector  are  sum¬ 
marized  in  Table  3,  along  with  the  equivalent  non-dimensionalized  parameter  set.  The  use  of  non- 
dimensional  parameters  using  the  Buckingham-71:  theorem  [84],  reduces  the  dimensionality  of  the 


FIGURE  27.  Snapshots  of  the  transient  injection  process  for  the  cross  injector  using  electrokinetic 
transport.  Starting  on  the  left,  the  analyte  in  red  is  drawn  from  the  sample  reservoir  or  upstream  mixer  at 
the  top  (1).  After  filling  the  injection  chamber  (2),  the  applied  potential  is  changed  to  inject  a  plug  of 
analyte  while  evacuating  leftover  sample  to  the  north  and  south  to  prevent  leakage  (3,4).  The  three  shorter 
channels  are  three  channel  lengths  long  and  the  long  channel  is  12  channel  lengths  long. 
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simulation  space  from  nine  independent  dimensions  (not  eleven,  since  the  values  for  the  electric 
fields  are  constrained  by  conservation  of  current)  to  four,  where  the  implicit  assumption  of  symmet¬ 
ric  shaping  (E2l  =  E. 4L >  Ej£>  =  E3D)  has  been  made. 

9  9 

The  analyte  fluidic  property  parameters  are  the  mobility,  //  (nr /Vs)  and  diffusivity,  /c(nr/s). 
The  stimuli  parameters  are  the  field  strengths,  EiL,  EiD  (V/m),  in  each  injector  leg  for  loading  (L) 
and  dispensing  ( D )  stages.  The  parameter  describing  the  injector's  topology  is  the  channel  width,  w 
(m),  with  an  implicit  assumption  that  all  channels  are  the  same  size.  The  non-dimensional  parame¬ 
ters  describing  the  operation  are  given  in  the  second  column;  they  are  the  ratio  of  accessory  fields  to 
driving  fields  in  loading  and  dispensing  stages,  ji1c,  n2c  and  the  Peclet  numbers  in  both  stages,  n3c, 

n4c- 

In  order  to  train  a  NN  model,  an  appropriate  set  of  model  outputs  must  also  be  selected.  In 
the  case  of  the  injector  device  the  output  is  the  set  of  parameters  that  describe  an  effective  Gaussian 
matching  the  variance  and  mass  of  the  plug  which  injected  into  the  separation  channel.  The  outputs 
shown  in  Table  3  are  the  output  peak  concentration  normalized  to  the  average  input  species  concen¬ 
tration,  and  the  output  band  variance  normalized  to  the  channel  width  squared. 

In  order  to  define  the  points  that  will  sample  the  design  space  defined  by  the  non-dimen¬ 
sional  input  parameters  of  Table  3,  we  first  sample  the  response  surface  along  traces  for  several 
nominal  values  of  the  input  parameters,  as  shown  in  Fig.  28a.  Several  useful  pieces  of  infonnation 
are  seen  in  these  graphs.  First,  the  Peclet  scale  is  exponential,  and  should  be  sampled  on  a  logarith¬ 
mic  scale.  The  field  ratio  scales  are  also  roughly  exponential  and  can  also  be  sampled  on  a  logarith¬ 
mic  scale.  The  scales  of  the  non-dimensional  parameters  are  chosen  to  encompass  a  large  region  of 
the  practically  realizable  injector  device  operation.  The  scales  are  listed  in  the  last  column  of 
Table  3.  With  the  information  of  the  design  space  traces,  the  entire  volume  can  be  sampled  using 
Gaussian  sampling.  The  results  of  this  sampling  for  various  projections,  two-parameters  at  a  time, 
are  shown  in  Fig.  28b.  The  lower  left  corner  of  the  7t4d  vs.  7t3d  plot  shows  points  that  were  moved 
out  of  the  infeasible  region  of  the  design  space.  The  total  of  900  simulations  were  used  to  sample 
the  design  space. 

In  order  to  select  the  appropriate  size  for  the  neural  network,  the  KFCV  technique  is  used  to 
measure  the  generalization  error.  A  plot  is  created,  shown  in  Fig.  29,  of  the  cross  validation  mea¬ 
sured  error  as  a  function  of  the  number  of  neurons  in  the  hidden  layer.  Since  the  Levenberg-Mar- 
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FIGURE  28.  (a)  Traces  through  the  design  space  for  fixed  parameters  of  ji1c  =  0.5,  n2c  =  0.5,  ji3c  =  400,  n4c 
=  400.  (b)  Scatter  plots  of  the  design  space  sampled  for  training  the  neural  network  model, 
quardt  algorithm  has  a  random  starting  point,  each  point  on  the  plot  is  the  average  of  1 5  repeated 

KFCV  tests.  The  calculations  for  the  cross  validation  and  neural  network  construction  are  done  in 

MATLAB.  Here  a  NN  is  created  with  the  topology  similar  to  that  in  Fig.  6b,  except  a  single  layer  of 

hidden  neurons  is  used.  The  hidden  layer  neurons  use  the  activation  function  of  Eq.  (10)  and  the 

output  neurons  use  the  linear  activation  function  in  Eq.  (11).  The  training  algorithm  is  the  Leven- 

berg-Marquardt  algorithm  [56],  for  a  maximum  of  350  iterations.  The  results  of  the  validation  test 

are  shown  in  Fig.  29.  The  result  shows  that  there  is  an  increasing  error  for  a  low  number  of  hidden 

nodes  due  to  large  biasing  error,  and  increasing  error  for  a  high  number  of  hidden  nodes  due  to 

larger  variance  errors,  with  a  minimum  near  50  nodes  for  both  models.  This  confirms  the  intuition 

that  simple  networks  will  have  trouble  fitting  all  of  the  original  data  points  due  to  lacking  degrees  of 

freedom,  whereas  highly  complex  networks  will  hit  the  original  data  points  easily,  but  vary  widely 

on  data  not  in  the  original  data  set,  due  to  too  many  degrees  of  freedom. The  noise  in  the  plot  of 

Fig.  29  decreases  with  an  increasing  number  of  averages  at  the  cost  of  additional  computation  time. 

Thus  this  type  of  plot  can  be  interpreted  to  say  that  on  average  networks  with  50  hidden  layer  neu- 
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FIGURE  29.  (a)  Mean  squared  error  results  from  cross  KFCV  model  selection  test.  The  minimum  occurs 
near  a  network  with  a  single  layer  of  50  neurons  for  both  models.  As  the  complexity  of  the  network 
increases,  the  generalization  error  increases  even  though  training  error  may  decrease  due  to  ‘overlitting’ 
the  data.  (b)The  trend  is  easier  to  see  with  the  data  logarithmically  scaled.  Noise  in  the  figures  is  from  the 
randomness  added  by  initializing  the  Levenberg-Marquardt  training  algorithm. 

rons  will  generate  a  model  with  the  best  generalization  error.  Thus  a  pool  of  networks  is  created  and 
trained  at  the  optimal  size,  and  the  network  with  the  lowest  training  error  for  our  model  is  selected. 

The  perfonnance  of  the  cross  injector  network  is  measured  by  comparing  to  a  validation  test 
set.  This  test  set  includes  simulations  in  the  7T3C-7r4c  plane  for  fixed  values  of  tt1c  =  0.9  and  n2c  = 
0.9.  Fig.  30  shows  the  results  of  the  neural  network  calculation  compared  to  numerical  solution  of 
the  governing  PDE's  to  determine  the  injected  band  variance  and  peak  height.  The  average  relative 
error  is  10.5%.  Increasing  the  pool  of  optimally  sized  networks,  or  adding  additional  data  points  to 
the  training  set  would  be  effective  ways  of  reducing  the  error  further. 

Double-Tee  Injector.  The  double-tee  injector  is  similar  to  the  cross  injector,  except  for  the  offset 
between  the  loading  channels.  The  most  common  value  for  this  separation  is  two  channel  widths, 
however,  for  this  example  this  value  is  a  variable.  The  results  of  an  injection  at  various  steps  in 
time,  for  two  channel  width  separation,  are  shown  in  Fig.  31  The  injector  operates  in  two  main 
stages,  loading  and  dispensing.  For  the  loading  stage  of  this  example,  the  analyte  shaded  in  red 
enters  the  injector  through  port  (A)  and  exits  through  port  (C)  into  a  waste  reservoir.  Simulta- 
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FIGURE  30.  Plots  of  band  variance  and  peak  height  vs.  n^c  for  varying  values  of  7i4c.  The  maximum  error 
of  both  models  is  approximately  10.5%.  Dots  are  the  numerical  simulations,  solid  lines  are  the  NN  function 
results. 


neously,  accessory  fields  on  the  sides  constrict  the  analyte  on  its  sides  from  ports  (B)  and  (D).  After 
the  analyte  has  completely  filled  the  injection  chamber,  the  applied  voltages  are  changed  for  the  dis¬ 
pensing  stage  so  the  electric  fields  push  the  analyte  plug  into  the  separation  channel  at  port  (B). 
While  the  main  driving  fields  push  the  plug  into  the  separation  channel,  accessory  electric  fields  are 
created  in  the  top  and  bottom  channels  at  ports  (A)  and  (C)  to  pull  analyte  back  from  the  injection 
intersection,  thereby  preventing  analyte  leakage  back  into  those  channels  while  the  plug  is  being 
injected. 

The  simulation  parameters  of  the  double-tee  model  are  found  in  Table  4  and  are  nearly  the 
same  as  the  cross  injector.  The  main  difference  is  the  double-tee  has  a  parameter  to  describe  the  off¬ 
set  of  the  input  channels  n This  parameter  is  varied  from  -2  to  2  to  represent  a  backward  to  for¬ 
ward  offset  of  two  channel  widths,  the  most  common  offset  found  in  the  double-tee  injector.  Fig.  32 
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FIGURE  31.  Snapshots  of  the  transient  injection  process  for  the  double-tee  injector  using  electrokinetic 
transport.  Starting  on  the  left,  the  analyte  in  red  is  drawn  from  the  sample  reservoir  at  the  top(l).  After 
filling  the  injection  chamber  (2),  the  applied  potential  is  changed  to  inject  a  plug  of  analyte  while 
evacuating  leftover  sample  to  the  north  and  south  to  prevent  leakage  (3,4,5). 
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Table  4. Double-Tee  Simulation  Parameters 


Simulation  Parameters 

System  Non-Dimensional 
Parameters 

Non-Dimensional  Parame¬ 
ter  Constraints 

Input 

Input 

Input 

/j  -  Electrokinetic  mobility 

7Ild  =  ^2L>/'£'lL 

rc3d,  n4d  e  [10,  5000] 

k  -Diffusion  coefficient 

"2d=  E]qJE2o 

nid’n2de  [0-01,5] 

T)l>  E2l,  E3l,  E4L-  Load¬ 
ing  electric  field 

"3d  =  /^3  1^/k 

n5de  [-2.2] 

Eu)’  AzD’  E4d  -  Dis¬ 

pensing  electric  field 

714(1=  MEi&w/ K 

3 

II 

tel 

w  -  Channel  width 

"5d=  L/w 

ElD~E3D 

C0  -  Analyte  concentration 

7t2d>  1/16.65  log(8490/ji4d) 

Output 

Output 

Cp  -  Equivalent  Gaussian 
Peak 

Cp  =  Cp/C0 

g2  -  Equivalent  Gaussian 
Variance 

g2  =  g2/w2 

shows  how  the  injector  topology  is  modified  with  this  parameter  and  its  affect  on  the  output  plug, 
all  other  parameters  held  constant.  As  the  parameter  becomes  equal  to  zero,  the  model  describes  the 
standard  cross  injector.  This  figure  shows  that  a  forward  offset  produces  the  largest  band,  followed 
by  a  zero  offset  and  then  the  backward  offset.  This  effect  is  due  to  the  manner  in  which  the  injection 
chamber  is  loaded,  and  the  magnitude  of  pullback  for  leakage  control.  Since  the  material  is  loaded 
from  the  top  (1)  and  pinched  from  the  sides  (2,4),  it  has  a  higher  concentration  towards  the  top  of 
the  channel,  as  seen  in  step  2  of  Fig.  31.  As  the  plug  travels  into  the  separation  channel  (2),  if  it 
must  pass  the  source  channel  (1),  it  will  lose  more  mass  than  if  it  must  pass  the  waste  channel  (3). 
This  is  due  to  the  vacuuming  effect  of  the  pullback  accessory  fields  in  channels  1  and  3,  and  the  fact 
that  there  is  more  mass  on  the  top  of  the  plug  than  on  the  bottom.  Thus  a  backward  offset  causes  the 


(a)  backward  offset  (b)  zero  offset  (c)  forward  offset 


FIGURE  32.  Simulations  of  the  double-tee  injector  for  jild  =  ji2d  =  2.5,  ji3d  =  n4d  =  500  for  (a)  Jtsd  =  -2  (b) 
jx5d  =  0  (c)  n5d  =  2.  The  surface  plots  are  not  to  scale,  (a)  and  (b)  are  scaled  to  a  maximum  concentration  of 
0.5  and  (c)  is  scaled  to  a  maximum  concentration  of  1.  (a)  is  a  backward  offset  which  produces  the  smallest 
plug,  (b)  is  a  zero  offset,  and  (c)  is  a  forward  offset  which  produces  the  largest  plug.  This  effect  is  apparent 
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FIGURE  33.  (a)  Trace  of  the  variance,  peak  height,  and  mass  vs.  channel  offset  for  7tj,j  =  =  0.5, 

n3d=400, 7t4d=400.  (b)  Trace  of  the  variance,  peak  height,  and  mass  vs.  channel  offset  for  7tld  =  7t2d  = 

2.5, 7t3d=400,  7t4d=400.  (c)  Projected  sampling  of  the  design  space. 

band  to  lose  more  mass  due  to  the  pullback  effect  and  a  forward  offset  causes  to  lose  the  least 
amount  of  mass  with  a  zero  offset  in  the  middle.  This  behavior  is  automatically  accounted  for  in  the 
numerical  simulations  used  to  train  a  neural  network  model,  and  is  less  significant  for  injections 
that  use  weak  pullback  fields. 

As  with  the  cross  injector,  traces  in  the  design  space  can  be  constructed  to  get  a  sense  of  the 
structure  of  the  response  surface.  The  operation  is  similar  to  that  of  the  cross,  except  for  the  offset 
parameter,  7r5(j.  At  nominal  values  of  7tlcj  =  7T2d  =  2.5,  7i3d  =  400,  n4(4  =  400,  the  plot  appears  in 
Fig.  33a.  When  the  effects  of  the  pullback  are  not  significant  7i2(j  <  1>  the  result  is  more  symmetric 
with  the  channel  offset,  as  seen  in  Fig.  33b.  Based  on  these  plots  7x5^  is  sampled  on  a  linear  scaling, 
and  all  other  parameters  are  sampled  on  a  logarithmic  scale,  as  with  the  cross  injector.  The  resulting 
sampled  space  is  seen  in  Fig.  33c.  A  total  of  4574  simulations  were  used  to  sample  the  space. 

Because  of  the  slightly  more  complicated  response  surface  structure,  a  neural  network  with 
two  hidden  layers  is  used,  as  in  Fig.  6b,  [44],  The  results  of  the  KFCV  model  selection  are  seen  in 
Fig.  34.  For  this  test,  the  two  layers  were  constrained  to  have  the  same  size,  leaving  one-degree  of 
freedom.  The  layers  could  be  sized  independently  to  create  a  two  dimensional  error  surface,  at  the 
cost  of  additional  computation.  The  increase  in  ‘overfitting’  error  is  again  evident  as  the  error 
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hidden  neurons 

FIGURE  34.  Mean  squared  error  results  from  double-tee  KFCV  model  selection  test.  The  minimum 
occurs  near  a  network  with  a  single  layer  of  35  neurons  for  the  variance  model  and  25  for  the 

decreases  and  then  increases  with  increasing  network  size.  The  optimal  network  size  for  the  vari¬ 
ance  model  occurs  near  35  nodes  in  each  layer  and  near  25  nodes  in  each  layer  for  the  peak  concen¬ 
tration  model. 

The  performance  of  the  double-tee  injector  network  is  measured  by  comparing  to  a  valida¬ 
tion  test  set.  This  test  set  includes  simulations  in  the  7t3d-7t4d  plane  for  fixed  values  of  7tld  =  0.9, 7t2d 
=  0.9,  and  7t5d  =  2.  Fig.  35  shows  the  results  of  the  neural  network  calculation  compared  to  numeri¬ 
cal  solution  of  the  governing  PDE's  to  determine  the  injected  band  variance  and  peak  height.  The 
maximum  relative  error  is  approximately  7.72%  for  both  the  variance  and  peak  concentration 
model.  Increasing  the  pool  of  optimally  sized  networks,  or  adding  additional  data  points  to  the  train¬ 
ing  set  would  be  the  most  effective  ways  of  reducing  the  error  further. 

Gated  Cross  Injector.  The  gated-cross,  Fig.  36,  has  the  same  geometry  as  the  cross,  but  uses  a 
three  stage  operating  scheme.  The  gated-cross  control  scheme  allows  for  continuous  flow  injection, 
so  that  new  sample  can  be  loaded  at  the  same  time  that  previously  dispensed  plugs  are  run  through 
the  separation  channel.  The  first  step  of  operation  involves  creating  the  gate  in  the  injection  cham¬ 
ber  by  counterflowing  an  analyte  stream  against  a  buffer  stream.  The  second  step  involves  remov¬ 
ing  the  applied  potential  from  the  buffer  port,  which  allows  a  portion  of  analyte  to  overflow  the 
injection  chamber.  The  final  step  is  a  return  to  the  applied  potentials  of  the  first  stage,  which  rees¬ 
tablishes  the  gate,  while  simultaneously  injecting  the  overflown  analyte  into  the  separation  channel. 
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FIGURE  35.  Plots  of  band  variance  and  peak  height  vs.  7t3c  for  varying  values  of  ji4c.  The  maximum  error 
for  both  models  is  approximately  7.72%.  Dots  are  the  numerical  simulations,  solid  lines  are  the  NN 
function  results. 

Since  the  action  taken  in  the  second  stage  involves  floating  only  one  node  while  all  others  remain 
unchanged,  no  independent  parameters  are  introduced  by  this  stage.  The  fields  and  flow  patterns  are 
completely  determined  by  the  parameters  set  in  the  first  stage. 

As  a  result  of  this  different  operation  the  gated-cross  has  a  set  of  parameters  that  differ  sig¬ 
nificantly  from  those  of  the  cross  and  double-tee,  as  seen  in  Table  5.  The  first  parameter,  rep¬ 
resents  the  extent  to  which  the  gate  is  closed.  As  discussed  in  [76],  as  long  as  Ej  >=E2,  the  gate  will 
remain  closed  in  the  limit  of  no  diffusion.  As  tt:1g  is  reduced,  the  gate  is  further  closed,  so  in  prac¬ 
tice  Ej  >  E2  in  the  presence  of  diffusion.  The  second  parameter  represents  the  ratio  of  the  buffer 


FIGURE  36.  Snapshots  of  the  transient  injection  process  for  the  gated-cross  injector  using  electrokinetic 
transport.  Starting  on  the  left,  the  analyte  in  red  is  drawn  from  the  sample  reservoir  at  the  top  and  takes  a 
90°o  turn  forming  a  gate  in  the  injection  chamber(l).  The  chip  operates  in  this  mode  until  a  plug  is 
needed.  At  this  time  the  electric  potential  in  channel  3  is  floated  causing  the  analyte  to  flood  the  chamber 
(2).  When  the  plug  of  desired  size  is  formed,  the  potential  is  reapplied  and  the  gated  is  reformed  injecting 
the  plug  into  the  separation  channel  (3). 
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electric  fields  (E3,  E4)  to  the  analyte  electric  fields  (E2,  E  | ).  As  k2q  increases,  the  buffer  fields 
become  larger  relative  to  the  analyte  fields.  The  third  parameter,  713(3  represents  the  Peclet  number 
of  the  loading  phase.  The  final  dynamic  parameter,  71:40,  is  a  measure  of  the  ratio  of  length  of  the 
injected  plug,  as  related  to  the  floating  time,  TLD,  to  the  channel  width. 

As  with  the  cross  and  double-tee,  a  significant  concern  in  the  operation  of  a  gated-cross 
injector  is  the  leakage  that  can  occur  in  the  separation  channel  if  the  gate  is  not  sufficiently  closed, 
as  seen  in  Fig.  38.  If  the  leakage  is  too  great,  the  injector  will  not  operate  because  the  increased 
noise  floor  will  make  a  separation  impossible.  A  region  of  feasibility  must  be  detennined  within 
which  the  injector  can  be  modeled.  This  leakage  was  analyzed  in  [76]  as  a  function  of  the  gate  clo¬ 
sure  and  the  system  Peclet  number.  They  detennined  a  boundary  indicated  by  1%  leakage  of  the 
flux  of  the  analyte  from  the  source  reservoir  into  the  separation  channel.  In  this  work,  we  determine 
the  boundary  for  a  more  complete  set  of  physical  parameters  to  create  a  region  of  feasibility. 

This  region  of  feasibility  is  determined  with  a  classification  neural  network.  For  classifica¬ 
tion  networks,  it  has  been  shown  that  a  decision  boundary  of  arbitrary  complexity  can  be  defined 
using  only  three  layers  (two  hidden,  one  output)  [44],  Using  the  simulations  like  those  in  Fig.  38, .a 
network  with  a  topology  of  Fig.  6b  with  two  hidden  layers  each  with  two  neurons,  is  created  to 
define  the  region  of  feasible  operation.Fig.  37  shows  the  feasibility  regions  for  the  gated-cross.  The 


Table  5. Gated  Cross  Simulation  Parameters 


Simulation  Parameters 

System  Non-Dimensional 

Non-Dimensional  Parame¬ 

Parameters 

ter  Constraints 

InDut 

InDut 

InDut 

/j  -  Electrokinetic  mobility 

nlg  =  E2,/El 

*13  e  [1/10,1] 

a- Diffusion  coefficient 

K2g=E3,/El 

JT2g  e  [L4] 

^il*  e2L’  ^3L>  e4L"  Load¬ 
ing  electric  field 

"3g  =  /jE3w/k 

nige  [10,  5000] 

w  -  Channel  width 

n4g=  w/juE2DTld 

n4gs  [1/20,  1/3] 

C0  -  Analyte  concentration 

OlltDUt 

OlltDUt 

Cp  -  Equivalent  Gaussian 
Peak 

Cp=Cp/C0 

a2  -  Equivalent  Gaussian 
Variance 

a2  =  a2/xv2 
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feasibility  space  is  independent  of  7t4G  because  the  leakage  is  determined  in  the  loading  stage  inde¬ 
pendently  of  how  the  band  is  dispensed.  The  most  notable  feature  of  the  feasibility  space  is  that  as 
713(3  the  Peclet  number,  becomes  smaller,  the  amount  of  feasible  space  decreases  significantly.  This 
is  due  to  the  fact  that  for  a  given  gate-closure,  when  diffusion  is  large  the  diffusive  flux  leaking  into 
the  separation  channel  is  large. 

The  gated-cross  model  is  sampled  in  a  similar  fashion  as  the  cross  and  double-tee.  The 
Peclet  number  is  sampled  on  a  logarithmic  scale,  all  others  are  linear.  The  points  in  the  infeasible 
region  are  moved  into  the  feasible  region  using  a  the  same  Gaussian  resampling  algorithm  with  the 
cross  and  double-tee.  For  a  network  with  a  single  hidden  layer  of  neurons,  the  resulting  optimal 
model  is  chosen  using  the  KFCV  algorithm,  Fig.  39a  The  optimal  model  sizes  are  approximately  35 
hidden  nodes  for  both  the  variance  and  peak  concentration  models.  As  a  qualitative  measure  of  the 


FIGURE  38.  Gated-cross  leakage  tested  at  Peclet  numbers  from  19  at  contour  1  to  oo  at  contour  5.  The 
contours  are  defined  by  the  7%  of  the  maximum  concentration,  and  agree  very  well  with  the  numerical 
and  experimental  results  found  in  [76] 
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accuracy  for  the  gated-cross  model,  the  response  surface  results  in  the  niG-n2Q  plane  are  shown  in 
Fig.  39b  and  Fig.  39c  for  the  variance  and  peak  concentration  models  respectively.  The  NN  model 
can  generate  a  much  denser  surface  of  points  in  a  fraction  of  second  compared  to  the  hours  it  takes 
a  numerical  simulation  to  generate  a  much  less  dense  surface.  This  type  of  computational  efficiency 
is  what  enables  simulation  and  synthesis  for  design,  where  many  repeated  simulations  may  be 
required. 

4.1.6.  Micro  Reactor  Model 

Traditional  chemical  manufacturing  is  heavily  based  on  economy  of  scale  with  large  reac¬ 
tors  and  associated  plants  requiring  large  process  batches  and  associated  large  scale  transport  and 
storage  of  raw  materials  and  products.  All  these  large  scale  features  present  health  and  safety  prob¬ 
lems  which  can  lead  to  major  disasters  as  well  as  unacceptable  levels  of  risk  to  operators  and  the 
neighboring  community  [85].  Microreactor  chemistry  shows  great  promise  as  a  novel  method  on 
which  to  build  new  chemical  technology  and  processes.  Micro  reaction  systems  are  fabricated  using 
techniques  originally  developed  for  electronic  circuits.  In  their  simplest  fonn,  micro  reactor  devices 
consist  of  network  of  micron-sized  channels,  typical  dimensions  are  in  the  range  of  30-300  pm, 


(a) 


ct2  model 


-3  C  model 

x  10  p 


(b) 

top  -  numerical  variance  model 
bottom  -  NN  variance  model 


(c) 

top  -  numerical  concentration  model 
bottom  -  NN  concentration  model 


FIGURE  39.  (a)  Mean  squared  error  results  from  gated-cross  KFCV  model  selection  test.  The  minimum 
occurs  near  a  network  with  a  single  layer  of  35  neurons  for  the  variance  and  concentration  models,  (b) 
Response  surface  of  numerical  simulation  on  top  and  neural  network  on  bottom  for  variance  model.  The 
neural  network  is  able  to  produce  results  in  a  fraction  of  section,  whereas  the  numerical  simulation  takes 
hours,  (c)  Response  surfaces  for  peak  concentration  model,  numerical  simulation  on  top,  neural  network 
on  bottom.  ji3q  and  ji4G  are  fixed  at  232  and  0.1208  respectively. 


56 


etched  into  a  solid  substrate  [85], [86], [87], [88],  For  solution-based  chemistry,  the  channel  networks 
are  connected  to  a  series  of  reservoirs  containing  chemical  reagents  and  products.  Reagents  can  be 
brought  together  in  a  specific  sequence,  mixed  and  allowed  to  react  for  a  specified  time  in  a  con¬ 
trolled  region  of  the  channel  network  using  electrokinetic  or  hydrodynamic  pumping.  We  consider 
electrokinetically  driven  systems  in  our  work.  In  these  systems  electrodes  are  placed  in  the  appro¬ 
priate  reservoirs  to  which  specific  voltage  sequences  can  be  delivered  under  automated  computer 
control.  This  control  offers  a  simple  but  effective  method  of  moving  and  separating  reactants  and 
products  within  a  microreactor,  without  the  need  for  moving  parts.  It  enables  the  ability  to  manipu¬ 
late  reagent  concentrations  in  both  space  and  time  within  the  channel  network  and  provides  an  addi¬ 
tional  level  of  reaction  control  which  is  not  attainable  in  bulk  stirred  reactors  where  concentrations 
are  generally  uniform.  Furthermore,  the  spatial  and  temporal  control  of  chemical  reactions  in  micro 
reactors,  coupled  with  the  features  of  very  small  reaction  volumes  and  high  surface  interactions,  is 
akin  to  the  situation  of  reactions  within  biological  cells. 

Much  of  the  current  research  on  microfluidic  reaction  systems  is  centered  on  the  fabrication 
of  micro  reactors  for  experimentation.  Ehrfeld  et.  al.  [86]  presented  an  overview  of  the  manufactur¬ 
ing  techniques  that  covered  the  fabrication  of  micro  mixers,  micro  reactors,  and  micro  heat 
exchangers.  Jensen  et.  al.  [89]  created  a  hydrodynamically  pumped  cross-flow  micro  reactor  and 
demonstrated  the  potential  of  it  as  a  laboratory  tool  for  heterogeneous  catalyst  testing.  Haswell  et. 
al.  [90]  considered  electrokinetic  based  fluidic  pumping  systems.  They  described  the  fundamental 
characteristics  and  emerging  applications  of  micro  reactors  in  the  field  of  synthetic  chemistry 
[85] [91],  While  the  micro  reaction  systems  are  undergoing  dramatic  expansions,  the  computer- 
aided  design  (CAD)  tools  have  not  kept  pace. 

In  this  project  we  developed  two  reduced  order  models  for  electrokinetic  continuous-flow 
micro  reactors.  The  models  can  take  any  concentration  profile  of  reactants  at  the  reactor  inlets  and 
simulate  diffusion,  convection,  and  reaction  inside  of  reaction  channels  of  any  length.  Both 
approaches  reduce  the  set  of  Partial  Differential  Equations  (PDEs)  to  a  system  of  coupled  Ordinary 
Differential  Equations  (ODEs).  Our  first  approach  is  based  on  physical  insights  and  consists  of 
approximating  the  reactive  channels  by  a  series  of  interconnected  parallel  Plug  Flow  Reactors 
(PFR).  The  second  approach  is  based  on  mathematical  simplification  of  the  underlying  equations. 
By  using  the  Method  of  Lines  (MOL),  the  spatial  derivatives  for  the  diffusional  mixing  in  reaction 
channels  are  approximated  by  finite  difference  relationships.  These  models  are  verified  by  numeri- 
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FIGURE  40.  Schematic  illustration  of  a  basic  micro  mixing  and  reaction  system. 


cal  simulations  in  COMSOL,  a  finite  element  solver  [92].  These  models  are  orders  of  magnitude 
faster  than  simulations  in  COMSOL,  and  therefore  allow  parametric  studies  of  a  particular  micro 
reactor  as  well  as  enable  integration  into  total  LoC  synthesis  and  layout  approaches  [93]. 

In  the  reaction  channel  as  shown  in  Fig.  40,  the  fluid  driven  by  electric  field  only  flows  in 
the  x  direction  and  the  species'  velocity  profiles  across  the  channel  are  flat  except  close  to  the  chan¬ 
nel  wall  [94]  [95].  As  a  result,  v,-  only  has  one  component  vjx  and  convection  happens  in  the  x  direc¬ 
tion  only.  Since  the  species  are  not  always  perfectly  mixed  before  entering  the  reaction  channel, 
their  concentrations  may  vary  in  both  the  x  direction  and  the  orthogonal  y  direction  across  the  chan¬ 
nel.  The  low  Reynold  number  operation  of  micro  scale  flows  [96]  indicates  that  mixing  in  the  y 
direction  is  caused  by  diffusion  only.  In  a  typically  designed  microfluidic  system,  the  fluid  flows 
fast  enough  that  the  concentration  changes  caused  by  diffusion  in  the  x  direction  is  much  less  than 
those  caused  by  convection.  In  channel  based  micro  reactors,  it  is  safe  to  only  consider  the  diffusion 
in  the  v  direction  and  ignore  the  diffusion  in  the  x  direction.  As  we  briefly  mentioned  earlier,  there 
are  two  operation  stages  for  LoC  system  and  mixing  and  reaction  are  belong  to  the  steady-state 
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FIGURE  41.  The  micro  reactor  is  represented  by  a  PFR  network  model.  The  calculation  of  reaction  and 
convection  is  conducted  in  one  column  of  PFRs  and  the  calculation  of  diffusional  mixing  is  conducted 
between  two  columns  of  PFRs. 


loading  stage.  We  further  assume  that  all  diffusion  coefficients  are  independent  of  concentration 
and  that  all  thermal  effects  are  negligible.  After  applying  these  assumptions,  we  have: 


2 

dC:  3  C; 

=  Di(—)  +  ri 

dx  dy2 


(38) 


Our  micro  reactor  models  are  based  on  this  set  of  PDEs.  Here  we  present  two  different 
approaches  that  approximate  this  PDE  system  as  an  ODE  system  which  is  easier  to  solve  and  nor¬ 
mally  has  shorter  computational  time. 


Plug  Flow  Reactor  Model.  In  our  microfluidic  system,  the  fluid  has  a  plug-like  velocity  profile. 
When  modeling  this  system,  it  is  natural  to  base  the  model  on  the  classic  Plug  Flow  Reactor  (PFR). 
To  capture  the  real  behavior  of  the  reactive  fluid  in  micro  reactors,  we  discretize  our  system  and  use 
multiple  PFRs  to  model  the  micro  reactor.  The  discretization  and  calculation  of  the  PFR  network 
model  is  shown  in  Fig.  41.  This  network  has  Mrows  and  N  columns.  Each  PFR  in  the  network  has 
the  length  of  /  and  the  width  of  w.  We  use  Pm  n  to  stand  for  the  small  PFR  in  the  /nth  row  and  nth 

column.  The  concentrations  of  the  species  inside  Pm  n  can  be  represented  as  cm  n‘(x),  where  the 
range  of  x  is  from  0  to  /.  The  total  volume  of  the  PFRs  in  the  network  remains  the  same  as  the  orig- 
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inal  micro  reactor  by  enforcing  =  w  and  =  L .  In  this  PFR  network  model,  the  calculation  of 

m  n 


concentration  change  is  conducted  alternatively  in  thex  and  y  direction. 

•  Step  0:  When  n  =  1,  calculate  the  initial  concentrations  for  the  first  column  of  PFRs  cm  /'(0)  by 
assuming  perfect  pre-mixing  for  each  PFR  and  average  the  inlet  concentrations  associated  with 

p 

1  m,l’ 

•  Step  1:  Eq.  (38)  with  the  diffusion  tenn  zeroed  out  is  used  to  calculate  convection  and  reaction 
inside  each  PFR  Pm  n  to  generate  cm  n‘(x).  When  n= 1,  the  initial  concentrations  cmJ{ 0)  are  cal¬ 
culated  by  Step  0;  otherwise,  they  are  set  to  the  values  of  the  modified  outlet  concentrations 
from  the  previous  column  of  PFRs  Pm>n.j : 


m,  n—  1 


(0 


•  Step  2:  Calculate  the  average  concentration  cm  n'  inside  Pm  n : 


(39) 


m,  n 


(40) 


This  average  concentration  is  used  in  the  next  step  for  approximating  the  diffusional  mixing. 


•  Step  3:  Fick’s  second  diffusion  law  is  used  in  the  y  direction  between  two  columns  of  PFRs  to 
approximate  the  concentration  change  caused  by  the  diffusional  mixing.  We  approximate  the 
second  derivative  with  finite  difference: 
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Vi 


(41) 


Step  4:  Calculate  the  modified  outlet  concentration  for  Pm  n : 


cm,n  (0  =  cm,n  (l)  +  Acm,n  (42) 

Then  the  calculation  will  repeat  from  Step  1  to  Step  4  until  reaching  the  final  column.  The 
combination  of  the  two  parameters  M  and  N  that  produces  accurate  results  for  a  system  with  given 


Peclet  number  ( Pe )  and  Damkohler  number  (Da)  is  not  unique.  A  fine  mesh  can  lead  to  a  result  with 
certain  accuracy,  but  requires  long  computational  time. 


Method  of  Lines  Model.  Eq.  (38)  is  the  governing  equation  for  our  micro  mixing-reaction  system. 
We  need  to  solve  the  set  of  PDEs.  In  this  section,  we  develop  the  micro  reactor  model  using  the 
numerical  method  of  lines  (MOL),  which  is  a  technique  for  solving  partial  differential  equations  by 
discretization  all  but  one  dimension,  and  then  integrating  the  semi-discrete  problem  as  a  system  of 
ODEs  [97].  A  significant  advantage  of  the  method  is  that  it  allows  the  solution  to  take  advantage  of 
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FIGURE  42.  Apply  the  numerical  method  of  lines  to  calculate  the  diffusion  in  the  channel  width  direction. 
Here,  we  take  10  discretized  concentration  points  in  the  channel:  c/  -  C;10.  Points  for  cf  and  c/1  are  needed 
to  apply  the  boundary  conditions. 


the  sophisticated  general  purpose  methods  and  software  that  have  been  developed  for  numerically 
integrating  ODEs.  It  is  necessary  that  the  PDE  problem  be  well-posed  as  an  initial  value  problem  in 
at  least  one  dimension.  In  our  micro  systems,  our  set  of  PDEs  satisfy  this  requirement. 

We  discretize  the  system  and  take  M  concentration  points  along  the  channel  width  as  shown 
in  Fig.  42,  in  which  M=10.  The  second  order  centered  fonnula  based  on  the  Taylor  expansion  can 
be  used  to  approximate  the  first  and  the  second  derivative: 

fix)  =  /(xA  +  1)~/(xA  1}  +  °ih2)  (43) 

2  h 

r(x)  =  /(**">-;;/'(/)+/ W1)  +  0(l2}  (44) 

n 

where  k  indicates  the  number  of  discretization  elements.  The  initial  conditions  for  the  set  of  ODE  are 
decided  by  assuming  perfect  pre-mixing  and  average  the  inlet  concentrations  for  each  discretization. 
As  the  micro  fluid  system  is  driven  electrokinetically,  the  boundary  conditions  for  the  fluid  flow  are 
treated  as  slip  boundary  conditions,  which  means  the  species  has  the  same  velocity  vix  at  each  point. 
The  Neumann  boundary  condition  [98]  is  used  for  the  concentration  change  at  the  channel  walls: 
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Table  6.  Timing  results  for  the  PFR  network  model. 


MXN 

2x2 

4x3 

4x4 

6x5 

6x9 

8x9 

10x9 

10x10 

error  (%) 

16.7 

8.8 

3.8 

2.6 

1.9 

0.7 

30.8 

0.2 

cpu  time  (s) 

0.058 

0.060 

0.072 

0.081 

0.107 

0.121 

0.143 

0.156 

-  1  %  M 

dct  ,  OC: 

—  =0  and  —  =0  (45) 

dy  dy 

By  using  Eq.  (43),  these  two  boundary  conditions  can  be  expressed  as 
=  c~  and  c^+l  =  cf  1 .  Then  the  original  set  of  PDEs  (Eq.  (38))  can  be  reformulated  into  a 
set  of  ODEs  as  follows: 
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(46) 


(47) 


(48) 


V  Ay'  2 

where  Ay  is  the  discretized  channel  width,  which  satisfies  =  w. 

m 

Increasing  the  amount  of  discretization  results  in  higher  accuracy  but  longer  computation 
time.  This  set  of  ODEs  can  readily  be  solved  by  a  standard  ODE  solver.  We  use  the  solver  called 
Lsode  to  solve  the  set  of  ODEs. 


Reaction  Model  Verification.  To  verify  the  accuracy  of  our  models,  we  developed  a  partial  differ¬ 
ential  equation  based  model  using  finite  element  methods  in  COMSOL  and  compared  the  results 
with  our  discretized  PFR  network  model  and  MOL  model.  The  case  study  involves  a  competitive 
reaction  system  with  species  A  and  B  competing  to  react  with  a  tag  T  to  produce  products  AT  and 
BT.  The  concentration  profiles  of  the  reaction  product  BT  generated  at  the  reactor  outlet  by  FEM 
simulation  is  shown  in  Fig.  43(a),  and  the  timing  results  as  a  function  of  element  size  are  shown  in 
Fig.  43(b). 

The  results  from  the  PFR  model  are  shown  in  Table  6.  It  takes  less  than  0.2  seconds  to  gen¬ 
erate  a  result  with  error  less  than  1%,  compared  to  the  few  minutes  required  for  similar  accuracy 
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FIGURE  43.  (a)  The  concentration  profiles  generated  by  FEM  simulations  with  decreasing  size  of  elements. 
The  red  curves  represent  the  results  associated  with  the  small  element  size.  They  are  grid  independent  and 
can  be  used  as  a  basis  for  comparison,  (b)  Timing  results  for  the  FEM  simulations. 


from  the  finite  element  simulation. 

The  results  for  the  method  of  lines  (MOL)  simulation  for  the  AT  product  concentration  and 
the  timing  results  are  shown  in  Fig.  44.  The  simulation  result  with  40  discretizations  takes  less  than 
0.25  seconds  to  generate  an  accurate  result,  i.e.  three  orders  of  magnitute  faster  than  the  FEM  simu¬ 
lation  (273  seconds). 


4.2.  System  Simulation 

This  section  demonstrates  the  integration  of  the  behavioral  models  for  subsystems  (e.g., 
mixers  in  Section  4.1.1  and  separation  microchips  in  Section  4.1.2.)  into  a  system-level  schematic 
of  integrated  LoCs.  Such  a  schematic  can  be  iteratively  simulated  to  capture  the  overall  effects  of 


Model  comparison.  FEM  simulation  vs.MOL  simulation.  Timing  results  for  MOL  model. 


FIGURE  44.  (a)  Comparison  of  MOL  model  simulation  results  with  FEM  simulations.  Using  more 
discretization  in  the  MOL  model  leds  to  a  more  accurate  solution,  (b)  Timing  results  for  the  MOL  model 
simulation  as  a  function  of  the  discretization. 
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subsystem-level  parameters  (e.g.,  separator  and  mixer  topologies  and  types;  operational  electric 
fields  and  injection  schemes  etc.)  and  element-level  parameters  (e.g.,  length  and  width  of  mixing 
and  separation  channels  etc.)  on  system  performance,  as  well  as  the  trade-offs  among  them,  leading 
to  a  system-level  optimal  design  achieved  by  “negotiation”  among  each  subsystem  and  element. 
Integrated  LoC  design  is  difficult  because  of  the  need  to, 

•  Accurately  interpret  the  fundamental  multi-physics  phenomena  (such  as  electric,  fluidic,  heat 
transfer  and  sample  transport)  at  element  or  component  levels. 

•  Properly  classify,  decouple  and  reduce  the  multi-physics  phenomena  to  tractable  mathematical 
models,  while  still  maintaining  their  essential  physical  characteristics. 

•  Efficiently  integrate  low-level  mathematical  models  to  obtain  a  system-level  representation  for 
iterative  simulation-based  design. 

Most  of  the  previous  approaches  to  LoC  modeling  are  not  amenable  to  accurate  simulation- 
based  iterative  design  of  integrated  systems  due  to  various  limitations,  such  as  difficulties  in  sys¬ 
tem-level  representation,  inaccurate  physical  descriptions  and  inflexibilities  with  respect  to  geomet¬ 
rical  perturbations.  Currently,  the  only  general  and  accurate  approach  is  the  classical  method  of 
numerical  simulation,  which  however  suffers  from  unacceptable  computational  cost.  Therefore,  a 
methodology  that  simultaneously  takes  into  account  of  balancing  accuracy  and  efficiency  without 
sacrificing  generality  is  needed. 

To  demonstrate  the  ability  of  the  methodology  developed  in  this  project  to  meet  these  goals, 
this  section  describes  how  it  can  be  used  to  design  a  competitive  immunoassay  microchip  consisting 
of  mixing,  reaction,  injection  and  separation  subsystems. 

4.2.1,  Integrated  Competitive  Immunoassay  Microchip 

Fig.  45a  illustrates  an  integrated  competitive  immunoassay  microchip  proposed  by  Harrison 

et  al.  [22],  which  consists  of  four  subsystems  (mixing,  reaction,  injection  and  separation).  Its  opera¬ 
tion  involves  both  synthesis  and  analysis,  which  are  typical  functions  performed  in  a  biochemical 
lab.  In  the  first  phase,  a  negative  electrical  voltage  is  applied  at  reservoir  3  with  reservoir  5,  6  and  7 
grounded  and  the  other  reservoirs  left  floating.  Theophylline  (Th,  from  reservoir  5)  and  fluorescein- 

labeled  theophylline  tracer  (Th  ,  from  6),  driven  by  the  electric  field,  first  mix  with  each  other 

within  channel  J1-J3.  Then  Th  and  Th  in  the  mixture  compete  for  a  limited  number  of  antibody 
(Ab,  from  reservoir  7)  binding  sites  in  the  mixing  and  reaction  channel  J3-J4.  The  mixture  (hereaf¬ 
ter  called  analyte)  of  reaction  products  Ab-Th*  and  Ab-Th,  as  well  as  unreacted  Th  and  Th  are 
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FIGURE  45.  (a)  Sketch  of  a  competitive  immunoassay  microchip  consisting  of  four  subsystems  (mixing, 
reaction,  injection  and  separation)  |22).  (b)  Its  system-level  schematic. 

loaded  into  the  double-T  injector.  Th,  Th  and  Ab  are  continuously  supplied  by  the  reservoirs; 
therefore  concentrations  of  all  the  samples,  reagents  and  products  in  the  mixer  and  reactor  at  this 
phase  are  in  steady-state.  This  completes  the  synthesis  operation. 

In  the  second  phase  (analysis),  the  voltages  are  switched  on  reservoirs  1  and  4  with  the  oth¬ 
ers  left  floating.  Thus  the  loaded  analyte  is  injected  as  a  band  into  the  separation  channel  14-reser¬ 
voir  4.  Because  the  analyte  includes  tagged  Ab-Th  and  Th  molecules  that  have  different  charges 
and  sizes,  they  move  at  different  speeds  and  eventually  can  be  separated  by  electrophoresis  [62]. 

The  amount  of  Ab-Th  and  Th*,  represented  by  the  areas  under  the  electropherogram,  for  a  wide 
range  of  Th  concentration,  can  be  obtained  to  generate  a  calibration  curve  to  determine  unknown 
concentrations  of  Th  in  the  sample  at  reservoir  5  for  clinical  analysis.  As  the  Ab-Th  complex  and 
Th  are  not  tagged  with  fluorescent  tracers,  they  are  invisible  to  detection  and  not  considered  in  the 

remaining  operation.  In  this  phase,  Ab-Th*  and  Th  bands  broaden  due  to  molecular  diffusion  and 
other  dispersion  sources;  therefore  the  transient  evolution  of  their  band  shapes,  directly  impacting 
the  analysis  quality  is  of  prime  importance. 

To  simulate  such  a  LoC,  a  system-level  schematic  that  assembles  behavioral  element  models 
based  on  the  LoC  geometric  and  functional  hierarchy  needs  to  be  built.  Fig.  45b  illustrates  the  hier- 
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FIGURE  46.  Behavioral  model  structure  for  the  electrokinetie  reactor.  At  the  inlet,  the  Fourier  coefficients 
( d of  widthwise  concentration  profiles  of  samples  and  reagents  conveyed  from  upstream  mixers 


characterize  their  premixing  degree.  At  the  outlet,  biofluidic  pins  quantify  the  average  concentrations  (c*™  *) 
of  the  reaction  products  and  unreaeted  samples. 


archical  decomposition  of  the  LoC.  It  is  decomposed  into  four  functional  subsystems  and  further 
broken  down  into  commonly  used  elements,  such  as  the  straight  mixing  channel,  injector,  reactor 
and  turn  separation  channel. 


4.2.2.  Behavioral  Models  for  Reactors  and  Injectors 

In  this  section,  behavioral  models  of  competitive  immunoassay  reactors  and  connection  pins 

(I/O  parameters)  of  double-T  injectors  that  are  integral  to  the  schematic  in  Fig.  45b  will  be  pre¬ 
sented. 

Competitive  Immunoassay  Reactors.  The  bio  fluidic  pins  for  a  general  purpose  reactor  will  be 
defined  to  enable  the  extension  of  our  system-level  simulation  approach  beyond  a  competitive 
immunoassay  design.  The  behavioral  model  for  the  reactor  that  will  be  developed  is  specific  to  the 
competitive  immunoassay. 

Pin  Definition.  In  practical  integrated  LoCs,  the  reactor  fulfills  two  functions,  bridging  the  mixer 
and  injector,  as  well  as  converting  samples  and  reagents  into  reaction  products  (synthesis).  There¬ 
fore,  the  biofluidic  pins  at  its  input  and  output  terminals  are  different.  Fig.  46  shows  the  behavioral 
model  structure  for  the  electrokinetie  reactor.  In  additional  to  electric  pins  (Vin  and  Vout),  bio  fluidic 
pins  are  also  proposed  with  arrows  indicating  the  direction  of  signal  flow  for  calculating  pin  values. 

At  the  reactor  inlet,  the  Fourier  coefficients  (d^)  of  the  widthwise  concentration  profiles  of  the 
samples  and  reagents  input  from  upstream  mixers  characterize  their  premixing  degree.  At  the  outlet, 

biofluidic  pins  are  defined  as  the  average  concentrations  (cj°“^)  of  the  reaction  products  and  unre¬ 
acted  samples  that  can  be  used  by  the  downstream  injector  to  determine  the  band  shape  of  the 
injected  species.  Here,  indices  in  and  out  represent  the  quantities  at  the  inlet  and  outlet  respectively. 
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FIGURE  47.  (a)  Sketch  of  the  competitive  immunoassay  reaction  between  Th,  Th  and  Ab.  (b)  Sketch  of  the 
competitive  immunoassay  reaction  model. 


Analog  wiring  buses  carrying  vector-type  pin  values  (concentration  coefficients  or  multiple  samples 
and  reagents)  similar  to  those  used  in  mixers  and  separators  are  also  employed. 

Behavioral  Model.  Fig.  47  depicts  the  process  of  competitive  immunoassay  reactions  between  Th, 
Th  and  Ab,  as  well  as  the  model  we  use.  The  uniformly  mixed  theophylline  Th  and  fluorescein- 

labeled  theophylline  Th  compete  for  a  limited  number  of  antibody  Ab  binding  sites  in  the  main 
mixing  and  reaction  channel  that  corresponds  to  channel  J3-J4  in  Fig.  45.  The  binding  kinetics  are 
governed  by  the  equations  in  Fig.  47a,  where  k\  and  kj  are  the  forward  and  backward  kinetic  con¬ 
stants  for  the  binding  reaction  between  Th  and  Ab  and  k 3  and  k4  are  those  between  Th  and  Ab 
respectively. 

The  concentrations  c  of  Th,  Th*,  Ab,  Ab-Th  and  Ab-Th  in  this  simultaneous  mixing  and 
binding-reaction  are  spatial-position  dependent  and  governed  by  a  set  of  convection-diffusion  equa¬ 
tions  with  reactive  source  terms  [99]: 
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where  D  is  the  sample  diffusivity,  u  is  the  sample  velocity  and  subscripts  Th,  Th*,  Ab,  Ab-Th  and  Ab- 
Th  refer  to  the  quantities  associated  with  the  corresponding  samples.  Eq.  (49)  is  non-linear  and  does 
not  admit  analytical  solutions.  Ref.  [22]  suggests  several  simplifying  assumptions: 


•  The  fluorescein  has  negligible  effects  on  the  binding  affinity  of  Th  to  Ab  (namely  kx  ~  k,  and 
k2  ~  k4),  and  on  the  material  properties  of  Th  and  Ab-Th  (namely  lln  ~  uTh-  and  Dn  ~  Dn>  ? 

U  Ab-Th  ~  U  Ab-Th"  aI1d  D Ab-Th  ~  &  Ab-Th"  )' 

•  The  reaction  is  essentially  irreversible,  k\  »  and  Ay  »  Ay,  such  that  the  forward  binding  is 
dominant  and  the  analyte  exiting  mixing-reaction  channel  does  not  contain  noticeable  Ab  (the 
amount  of  Ab  is  limited).  This  can  be  inferred  from  the  electropherogram  (Figure  6)  in  Ref.  [22] 

(if  the  reverse  reaction  was  significant,  another  concentration  peak  of  Th  disassociated  from 

Ab-Th*  complex  would  be  observed  in  the  separation  channel). 

•  Both  mixing  and  reaction  are  complete,  which  are  implied  by  Figures  5  and  7  in  Ref.  [22]  (Fig¬ 
ure  5  shows  that  99  %  mixing  has  been  reached.  In  Figure  7,  except  at  Th  concentration  of  10 
mg/L,  the  stop-flow  case  does  not  appreciably  enhance  the  reaction,  indicating  that  the  reaction 
is  almost  complete  for  the  continuous-flow  case). 

Based  on  these  assumptions,  Eq.  (49)  can  be  reduced  to 
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The  first  two  equations  in  Eq.  (50)  differ  from  each  other  by  a  scaling  constant  CR ,  the  con¬ 
centration  ratio  of  cT/r  to  cj y7*,  implying  that  the  amount  of  Ab  bound  to  Th  and  Th  (namely,  the 


c Ab-Th  anc'  c Ab-Th*)  is  proportional  to  the  concentrations  of  Th  and  Th  and  CR  holds  unchanged 
from  the  inlet  to  the  outlet  of  the  mixing-reaction  channel.  Additionally,  the  assumption  of  complete 
mixing  and  binding-reaction  enables  the  calculation  of  the  concentrations  of  reaction  products  and 
unreacted  samples  at  the  end  of  the  mixing-reaction  channel  (intersection  J4)  based  on  the  mass- 
balance  principle.  In  the  model  implementation,  the  main  mixing-reaction  channel  in  Fig.  47a  is 
modeled  as  a  serial  connection  of  a  mixing  channel  having  the  same  dimensions  as  the  original 
channel  and  a  reactor  with  negligible  physical  size  (R  =  0)  in  which  the  binding  reaction  completes 
instantaneously  (Fig.  47b).  This  treatment  allows  the  designer  to  monitor  the  mixing  degree  (or 
mixing  residual)  before  samples  and  reagents  enter  the  reactor  model  and  examine  if  the  complete 
mixing  assumption  is  satisfied  during  design  optimizations. 

The  mixing  model  is  described  in  the  previous  chapter  and  will  not  be  repeated  here.  The 
modeling  effort  hence  focuses  on  the  reactor  in  Fig.  47b.  As  the  amount  of  Ab  is  limited  (eventually 
completely  consumed),  the  binding  reaction  in  the  reactor  model  can  be  quantitatively  described  by, 
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Table  7. Molar  mass  flow  rates  of  samples  and  reaction  products  at  the  inlet  and  outlet  of  the  reaction 
model 
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where  =  uA  ■  c1"'1  • A  •>  Ar  is  the  constant  cross-sectional  area  of  the  mixing-reaction  channel,  M 
stands  for  the  molar  mass  flow  rates  of  samples  and  indices  in  and  out  represent  the  quantities  at  the 
inlet  and  outlet  of  the  reactor  model.  Similarly,  those  of  Th  and  Th  are  expressed  as 
=  u.r,  .c(i:]  ■ A  ,  M{m)  =u  ,  ■  c{m)  ■ A  .  The  contribution  of  the  axial  diffusive  mass  flux  is  not 
counted  due  to  the  long  mixing-reaction  channel  (L  =  81  mm  and  w  =  52  pm).  Then,  from  mass  bal¬ 
ance,  the  molar  mass  flow  rates  of  the  reaction  product  Ab-Th  and  unreacted  Th  exiting  the  reactor 
model  are  respectively  given  by, 
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Table  7  summarizes  the  molar  mass  flow  rates  of  samples  and  reaction  products  flowing  in 
and  out  of  the  reactor. 
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FIGURE  48.  Pins  for  the  electrokinetie  injector.  Indices  in,  out,  b,  aw  refer  to  the  quantities  at  the  terminals 
respectively  linking  to  the  reactor,  separation  channel,  buffer  reservoir  and  sample  (analyte)-waste  reservoir. 


where  uA^  ~  UAh-Th*  is  used  [22].  The  uniform  species  concentrations  (namely  the  average  concen¬ 
tration  attributed  to  complete  mixing)  ,  c^')  and  4“  '  at  the  inlet  of  the  reactor  model  can  be 
extracted  from  the  zeroth  component  of  .  Eqs.  (55)  and  (56)  establish  the  signal  flow  relationship 
between  the  species  concentrations  at  the  inlet  and  outlet  of  the  reactor  model. 

Obviously,  this  model  does  not  provide  reaction  kinetics  and  requires  sufficiently  long  chan¬ 
nels  to  ensure  both  complete  mixing  and  reaction  (this  assumption  is  valid  for  the  integrated  LoC  in 
Fig.  45),  leading  to  conservative  assessment  of  the  mixing-reaction  channel  length  for  LoC  design. 

Double-T  Injectors.  The  double-T  injector  [22]  serves  as  a  physical  junction  between  the  reactor, 
separator,  sample  (analyte)-waste  and  buffer  feed  channels.  It  operates  at  both  loading  and  dispens¬ 
ing  phases  and  introduces  analytes  from  the  continuous-flow  reactor  into  the  separation  channel  that 
involves  transient  evolution  of  species  bands.  Therefore,  it  is  indispensable  in  transitioning  the  LoC 
from  the  synthesis  to  analysis  phases.  In  contrast  to  the  injector  model  used  for  separation  modeling 
(Section  4.1.2.),  this  new  model  is  capable  of  automatically  detennining  the  shape  of  the  injected 
species  band  based  on  the  electric  fields  at  both  phases  and  analyte  concentrations  from  the  adjacent 
reactor  model. 

Pin  Definition.  Fig.  48  shows  the  pins  at  the  injector  terminals.  Since  the  injector  is  modeled  as  a 
single  element,  the  symbol  view  does  not  reflect  its  real  physical  geometry,  e.g.,  double-T  or  cross 
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[100].  The  terminals  connecting  to  the  reactor  and  separation  channel  are  respectively  defined  as 
the  inlet  and  outlet,  indicating  the  direction  of  signal  flow  for  calculating  the  band  shape  of  the 

injected  species.  At  the  inlet,  the  information  of  the  average  species  concentrations  (cj^)  is 
acquired  from  the  upstream  reactor  model.  At  the  outlet,  the  initial  species  band  shape,  such  as 
skew  coefficients  (S[ou>)),  variance  (a2  )  and  amplitude  (Aout),  as  well  as  the  starting  separation 
time  (tout  =  0)  is  set  and  propagated  to  the  downstream  separation  channel  for  separation  computa¬ 
tion.  In  addition,  four  electric  pins  Vin,  Vout,  Vh  and  Vaw  are  also  defined,  where  indices  in,  out,  aw 
and  b  represent  the  quantities  at  the  terminals  linking  respectively  to  the  reactor,  separation  channel, 
sample  (analyte)-waste  feed  channel  and  buffer  feed  channel. 

Behavioral  Model.  As  its  flow  path  lengths  are  negligibly  small  compared  with  those  of  feed,  mix¬ 
ing-reaction  and  separation  channels,  injector  can  be  assumed  as  a  point- wise  entity  and  electrically 
represented  as  four  resistors  with  zero  resistance  between  each  tenninal  and  the  internal  node  A), 

Kw  =  Rb  =  Rin  =  Kut=  o  (57) 

Here,  Nj  corresponds  to  the  intersection  of  flow  paths.  Thus,  the  voltages  at  the  terminals  are 

consequently  the  same  (y  =  y  =  y  =V  )•  The  functional  relationship  of  the  biofluidic  states 
between  the  inlet  and  outlet  of  the  injector  is  detennined  by  electric  fields  in  the  adjacent  channels 
and  quantitatively  described  by  a  parameterized  reduced-order  model  developed  with  the  neural 
network  approach  [40]. 

4.2.3.  Connection  of  Pins 

With  all  behavioral  models  available,  the  next  step  is  to  link  pins  at  element  terminals  by 
wires  (electric)  and  wiring  buses  (biofluidic)  according  to  the  spatial  chip  layout  and  compatibility 
of  the  physics  to  accomplish  the  simulation  schematic. 

In  Table  8,  the  bio  fluidic  wiring  buses  are  classified  into  the  intra-subsystem  (connecting  the 
element  within  a  subsystem,  e.g.,  mixing  or  separation)  and  inter-subsystem  (linking  the  sub¬ 
systems  involving  different  functions)  buses.  Since  the  reactor  and  injector  are  modeled  in  tenns  of 
a  single  element,  there  is  no  intra-subsystem  pin  definition. 
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Table  8. Biofluidic  wiring  buses  used  in  the  integrated  LoC  simulations 


Connection  Type 

Subsystem  Type 

Bus 

Pin  Name 

Mixing 

d  [0:291 

Concentration  coefficients 

Intra-subsystem 

Separation 

t  [0:2] 

o2  [0:2] 

A  [0:2] 

5  [0:  301 

Separation  times 

Variance 

Amplitude 

Skew  coefficients 

Reaction 

N/A 

N/A 

Injection 

N/A 

N/A 

Mixing-Reaction 

d  [0:291 

Concentration  coefficients 

Reaction-Inj  ection 

cflVg[0:2] 

Average  analyte 
concentrations 

Inter-subsystem 

Injection-Separation 

t  [0:2] 
a 2  [0:2] 

A  [0:2] 

5  TO:  301 

Separation  times 

Variance 

Amplitude 

Skew  coefficients 

4.2.4.  System  Simulation  of  Integrated  LoC  Systems 

In  this  section,  the  system-level  schematic  simulation  procedure  will  be  first  described. 

Then,  the  simulation  results  of  the  integrated  competitive  immunoassay  microchip  will  be  discussed 
and  validated  against  numerical  and  experimental  data. 

Simulation  Description.  The  immunoassay  schematic  is  simulated  in  two  consecutive  steps  arising 
from  the  two  operational  phases.  For  each  step,  both  electric  and  biofluidic  simulations  are  con¬ 
ducted. 

•  Mixing-reaction-loading  (synthesis)  phase.  In  this  phase,  voltage  is  applied  at  reservoir  3  with 
reservoirs  5,  6  and  7  grounded  and  the  others  left  floating  (potential  setting  is  inactivated  in  the 
schematic).  Given  system  topology  and  element  dimensions,  nodal  voltages  at  element  termi¬ 
nals  within  the  entire  system  are  first  computed  by  Ohm’s  and  Kirchhoflf’s  laws  using  the  resis¬ 
tor  models.  The  resulting  nodal  voltage  and  current  through  the  element  are  in  turn  used  to 
calculate  the  electric  field  strength  ( E )  and  its  direction  within  the  mixer,  reactor  and  analyte- 
waste  feed  channel,  as  well  as  the  flow  ratios  at  intersections  J1  and  J3.  With  these  results  and 
user-input  sample  properties  (D  and  //),  the  biofluidic  pin  values  are  computed  sequentially 
starting  from  the  most  upstream  sample  reservoirs  (5,  6  and  7)  to  the  injector.  The  sample  (Th, 
Th  and  Ab)  concentrations  after  mixing  are  determined  and  then  fed  to  the  reactor  model  to 
calculate  the  concentrations  of  detectable  species,  e.g.,  the  reaction  product  (Th  -Ab)  and  the 

unreacted  Th*.  At  the  double-T  injector,  the  concentration  infonnation  along  with  the  electric 
fields  in  the  adjacent  channels  is  saved. 

•  Dispensing-separation-detection  (analysis)  phase.  In  this  phase,  the  voltage  is  switched  to  reser¬ 
voirs  1  and  4  with  the  others  left  floating.  As  the  synthesis  phase,  the  nodal  voltage  in  the  net¬ 
work  is  first  computed.  Then,  the  computation  of  separation  state  (e.g.,  the  arrival  time, 
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Table  9.Comparison  between  system-level  simulation  results  and  experimental  data  on  normalized  concentrations  of 
Th+  at  channels  J1-J3  and  J3-J4. 


Channel 

Svstem  Simulation 

Experiment 

J1-J3 

0.5 

0.534 

J3-J4 

0.25 

0.257 

variance,  skew  and  amplitude)  is  initiated  at  the  injector  outlet  using  the  stored  infonnation 
from  phase  (1),  and  subsequently  propagated  to  downstream  separation  channels  using  signal 
flow  until  the  waste  reservoir  is  reached. 

Simulation  Results  and  Discussion.  In  this  section,  the  system  simulation  results  and  their  com¬ 
parison  to  experimental  and  numerical  data  will  be  presented  to  validate  the  modeling  methodology. 
Table  9  shows  the  comparison  between  the  system  simulation  results  and  experimental  data 

on  the  normalized  concentrations  of  Th  (by  the  initial  reservoir  sample  concentration)  at  channels 
J1-J3  and  J3-J4.  Since  channels  J1 -reservoir  5  and  J1 -reservoir  6  have  exactly  same  dimensions 
(that  is  same  electric  resistance),  flow  rates  through  them  are  also  identical  with  a  single  constant 

potential  applied  at  reservoirs  5  and  6.  This  leads  to  a  two-fold  dilution  of  Th*  at  channel  J1-J3. 
Likewise,  channel  J3-J2-reservoir  7  is  fabricated  with  the  same  resistance  as  the  combined 
resistance  of  channels  J3-J1 -reservoir  5  and  J3-J1 -reservoir  6,  eventually  yielding  a  four- fold 
dilution  of  Th*  at  channel  J3-J4.  The  simulation  results  agree  with  the  experimental  data  excellently 
with  the  worst  error  of  7  %. 

In  Fig.  49,  calibration  curves  of  the  area  ratio  for  unreacted  Th*  and  Ab-Th*  complex  are 
obtained  from  system  simulation  results  by  varying  the  concentration  of  Th  (Th*  concentration  is 

fixed).  Area  ratios  of  Th  and  Ab-Th  are  respectively  defined  as  Ar^, +  Ar^  n;  )  and 

Ar ih_n*  j ( Ar/(,,  +  Ar^  n, ) ,  where  Ar  is  the  species  amount,  represented  by  its  area  under  the  elec- 
tropherogram  (Fig.  50).  The  simulation  results  match  the  experimental  data  very  well  with  a  rela¬ 
tive  error  less  than  5  %.  As  the  concentration  of  Th  increases,  the  amount  of  Th  begins  to 

predominate  in  the  mixture  of  Th  and  Th  in  mixing-reaction  channel  J3-J4.  Therefore,  Ab  exhibits 
a  higher  probability  of  colliding  and  binding  to  Th  than  to  Th*,  leading  to  more  unreacted  Th  and 
the  growth  of  Th*  area  ratio. 

Fig.  50  shows  the  simulated  electropherograms  of  separating  Th*  and  Ab-Th  at  three  detec¬ 
tion  spots  in  the  separation  subsystem.  It  is  seen  that  species  bands  of  Th  and  Ab-Th  gradually 
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FIGURE  49.  Comparison  between  the  system-level  simulation  results  and  experimental  data  (22]  on  the  area 
ratio  of  unreaeted  Th*  and  Ab-Th*  complex.  Abscissa  shows  the  initial  concentration  of  Th  before  the  50- 
fold  off-chip  dilution.  Actual  concentrations  of  Th*  and  Ab  are  not  available.  Their  values  were  extracted 
from  the  results  in  Ref.  |22]. 


Table  10. Comparison  between  numerical  analysis  [29]  and  system  simulations  on  the  area  ratio 


Numerical 

Sys.  Simul 

Area  ratio  Th 

0.845 

0.84 

Area  ratio  Ab-Th 

0.155 

0.16 

separate  during  their  migration  through  the  separation  channel  but  with  a  considerable  decrease  in 
band’s  amplitude  and  spreading  in  band’s  width  due  to  the  dispersion.  The  area  ratio  of  the  species 
bands  are  extracted  at  the  third  detector  (bottom  trace)  and  compared  with  numerical  simulation 
[29]  in  Table  10.  Very  good  agreement  with  the  worst-case  error  less  than  10  %  and  impressive 
speedup  >100x  over  the  numerical  simulations  have  been  achieved. 

4.2.5.  An  Improved  Design 

Finally,  the  system-level  simulation  is  employed  to  redesign  the  original  chip  and  improve 
its  bio-analysis  efficiency  and  minimize  its  chip-area.  Here,  the  original  constraints  on  the  channel 
size  and  power  supply  are  kept.  Specifically,  the  mixing,  reaction,  injection  and  separation  channels 
are  fabricated  with  channel  width  52  pm  and  those  feed  channels  leading  to  buffer  and  sample  (ana- 
lyte)-waste  reservoirs  are  236  pm  wide.  The  same  power  supply  with  a  single  output  of  6  kV  is 
used.  Due  to  the  practical  fabrication  limit,  I/O  reservoirs  are  required  to  be  spaced  at  least  5  mm 
apart.  In  addition,  the  immunoassay  reaction  is  also  assumed  to  be  mixing-limited,  thus  complete 
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FIGURE  50.  Electropherograms  of  unreacted  Th  and  Ab-Th  complex  from  system  simulations  at  three 
detection  spots:  10  mm  after  injection  (top  trace),  before  the  90°  elbow  (middle  trace)  and  after  the  90°  elbow 
(bottom  trace).  Th  concentration  used  in  experiments  and  simulations  is  40  mg/L  before  the  50-fold  off-chip 
dilution  (800  pg/L  after  the  dilution). 

mixing  also  signifies  complete  reaction  (otherwise  a  period  of  incubation  needs  to  be  added  at  the 
end  of  the  synthesis  phase  from  the  design,  in  which  all  reservoirs  are  left  floating). 

Fig.  51  shows  the  redesign  by  reducing  the  excessive  mixing  length  and  arranging  it  into  a 
compact  serpentine  geometry.  The  original  mixing  channel  J1-J3  and  J3-J4  are  shrunk  from  26.5 
and  81.63  mm  to  9.3  and  77.33  mm  respectively,  while  keeping  almost  the  same  mixing  degree. 
Additionally,  the  overly  long  separation  channel  is  shortened  from  81.11  to  25.5  mm,  which  hence 
increases  the  electric  field  (from  600  V/cm  to  1.8  kV/cm)  and  minimizes  the  band  spreading  due  to 
diffusion  (Joule  heating  dispersion  is  still  negligible  in  this  case).  Furthermore,  the  detector  is 
moved  to  the  front  of  the  90°  elbow  separation  channel  to  avoid  the  turn  dispersion  at  the  high  elec¬ 
tric  field  in  the  new  design. 

Table  11  shows  the  modified  mixing  and  separation  channel  length,  as  well  as  the  system 
performance  between  the  original  and  improved  designs.  Although  the  separation  resolution  drops 
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FIGURE  51.  Schematic  (not  to  scale)  of  the  improved  design  of  Fig.  45  using  the  system-level  simulation, 
to  15,  it  is  still  high  enough  to  resolve  the  species  bands.  Most  importantly,  an  impressive  1.6-fold 
increase  in  concentration  amplitude  and  nearly  3-fold  and  10-fold  decreases  in  variance  and  chip 
area  have  been  achieved.  In  other  words,  an  assay  chip  that  is  easier  to  detect  has  been  designed  in 
less  area. 

4.2.6.  Summary 

This  chapter  describes  the  system-level  schematic  simulation  of  an  electrokinetic  competi¬ 
tive  immunoassay  LoC.  In  contrast  to  the  separation  and  mixing  models  from  the  previous  chapters, 


Table  1 1  .Comparison  of  the  channel  dimension,  separation  time,  variance,  peak  height  and  chip  area  between  the 
original  design  and  the  improved  design. 


Parameters 

Original 

Improved 

J1-J3 

26.5  mm 

9.3  mm 

J3-J4 

81.6  mm 

77.33  mm 

Mixing  degree 

99% 

99% 

J4-reservoir  4 

81.11  mm 

25.5  mm 

Ab-Th* 

Th* 

Ab-Th* 

Th* 

Separation  Time 

18.6 

29.1 

2.12 

3.32 

(s) 

Variance  (pm2) 

66607 

30091 

15926 

11400 

Amplitude  (norm.) 

1 

4.9 

2.4 

8 

Resolution 

24 

15 

Chip  area 

7.6  cm  x  7.6  cm 

2.5  cm  x  2  cm 
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pins  involving  different  physics  have  been  defined  at  the  tenninals  of  reactor  and  injector  models, 
which  enable  the  interconnection  among  different  subsystems  and  transition  between  the  opera¬ 
tional  phases  of  the  chip.  A  simplified  reactor  model  valid  for  forward  immunoassay  reactions  has 
been  developed,  implemented  in  an  analog  hardware  description  language  (Verilog-A)  and  linked  to 
behavioral  models  of  mixers,  separation  channels  and  injectors  [40]  to  fonn  a  complete  simulation 
schematic.  An  overview  of  the  pin  connections  within  and  between  subsystems  has  been  provided 
to  interpret  the  biofluidic  signal  transmission  within  the  entire  LoC  network. 

The  constructed  schematic  then  has  been  used  to  simulate  the  LoC  in  both  synthesis  and 
analysis  phases,  each  requiring  sequential  solutions  of  the  electric  network  by  Kirchhoff’s  law  and 
biofluidic  states  by  signal  flow.  The  simulation  results  of  on-chip  mixing,  reaction,  injection  and 
separation  have  been  verified  by  numerical  and  experimental  data.  A  speedup  (>100x)  is  achieved 
over  the  numerical  FEM  simulations,  while  still  maintaining  high  accuracy  (relative  error  less  than 
10  %).  The  system-level  simulation  captures  the  overall  effects  of  chip  topology,  element  size  and 
operational  parameters  on  chip  performance  and  is  used  to  redesign  the  original  chip  to  improve 
analysis  quality  but  occupy  less  chip-area.  This  application  of  the  modeling  and  simulation  frame¬ 
work  developed  in  this  thesis  demonstrates  its  effectiveness  in  system-level  design  of  integrated 
LoCs. 

4.3.  Component  Optimization 

Our  design  approach  has  three  fundamental  components:  (1)  The  system-level  simulator 
described  in  the  preceding  sub-sections,  which  is  capable  of  simulating  complex  designs  in  seconds 
(2)  General  area  constraints  and  heuristic  Channel  Packing  and  Placement  algorithms  (CPP).  (3) 
Tailored  design  optimization  methods  that  utilize  the  system  simulator,  CPP  algorithms,  numerical 
optimization  and  area  constraints. 

4.3.1.  Design  Constraints 

We  have  approached  the  design  of  micro  electrophoretic  devices  from  two  perspectives. 
Specifically  we  wish  to:  (1)  Determine  the  highest  perfonnance  design  that  can  fit  within  a  given 
chip  area.  (2)  Minimize  the  area  occupied  by  a  design  while  maintaining  device  feasibility.  The 
trade-offs  between  these  two  objectives  will  be  discussed  in  this  report. 

We  have  divided  the  constraints  on  the  design  problem  into  performance  constraints,  oper¬ 
ating  constraints  and  physical  constraints  as  summarized  in  Table  12. 
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Table  12. Examples  of  three  different  constraint  classes  along  with  typical  numerical  values. 


Perfonnance 

Operating 

Physical 

N  >  104 

C>  lpM 

co  e  [10,100] 

pm 

H  <  10  mm 

E  <  10  V/cm 

Area  <  25  cm2 

R>  1.5 

V  <  30  kV 

PAD  >  5  pm 

Performance  constraints  provide  a  lower  bound  on  the  effectiveness  of  the  separation  sys¬ 
tem.  In  our  case  we  require  that  separations  always  achieve  baseline  resolution  (>  1.5).  Operational 
constraints  are  dictated  by  external  equipment  such  as  the  available  voltage  source  or  the  detector 
sensitivity.  Even  if  this  equipment  does  not  occupy  physical  space  on  a  chip,  it  can  have  a  signifi¬ 
cant  impact  on  the  final  design.  Physical  constraints  are  those  dictated  by  regions  of  model  applica¬ 
bility  and  fabrication  limits.  Fabrication  constraints  capture  manufacturing  limits;  for  instance,  the 
minimum  channel  width  co  or  the  minimum  spacing  between  channels  PAD. 

Currently,  since  the  range  of  possible  channel  widths  co  and  channel  depths  5  is  small,  we 
treat  these  quantities  as  parameters  and  fix  their  values.  In  general,  we  find  that  smaller  channel 
widths  provide  better  perfonnance  and  smaller  area.  Furthermore,  our  models  can  neglect  the 
boundary  layer  effects  that  become  important  in  very  small  channels.  Our  smallest  channel  dimen¬ 
sion  is  well  over  an  order  of  magnitude  larger  than  the  Debye  length  [36]  produced  by  a  standard 
CE  electrolyte  in  a  glass  or  plastic  chip. 

Connectivity  and  Area  Constraints.  The  two  most  important  physical  constraints  are  what  we  call 
connectivity  constraints  and  area  constraints.  Connectivity  constraints  are  used  to  uniquely  specify 
a  particular  channel  topology.  Area  constraints  are  used  to  detennine  the  space  requirement  for  a 
given  channel  topology  on  the  chip. 

Fig.  52  illustrates  the  three  parts  of  our  connectivity  constraints  by  showing  the  infonnation 
necessary  to  connect  a  U-bend  to  a  straight  channel  section.  It  is  apparent  that  the  U-bend  must 
attach  to  the  end  of  the  straight.  However,  this  results  in  four  possible  scenarios  as  shown  in 
Fig.  52(a).  Two  of  these  turns  are  illegal  (turns  3  and  4)  because  they  double  back  on  the  straight 
section.  Doubling  back  can  be  eliminated  by  enforcing  a  flow  direction  as  shown  in  Fig.  52(b).  The 
flow  direction  is  a  cardinal  direction  {north,  south,  east,  west}  where  the  flow  direction  out  of  an 
upstream  section  must  equal  the  flow  direction  in  to  the  downstream  section.  This  leaves  turns  1 
and  2  as  viable  alternatives.  The  topology  can  be  uniquely  specified  by  providing  a  rotational  direc- 
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FIGURE  52.  The  three  necessary  connectivity  constraints,  (a)  Two  legal  turns  (1  and  2)  and  two  illegal  turns 
(3  and  4)  that  can  be  attached  to  the  endpoint  of  a  straight  channel,  (b)  The  outlet  flow  direction  of  the 
straight  and  the  inlet  flow  direction  of  the  turn  must  match,  thereby  eliminating  illegal  turns,  (c)  The 
topology  becomes  unique  when  a  turn  direction  is  specified. 

tion  {clockwise,  counterclockwise}  to  the  turn  as  shown  in  Fig.  52(c).  These  three  pieces  of  infor¬ 
mation  allow  a  designer  to  build  a  unique  2  section  channel  topology  from  the  library  of  section 
types  discussed  in  Section  4.2. 

However,  the  amount  of  information  required  to  uniquely  specify  any  particular  topology 
grows  with  the  number  of  channel  sections.  Fig.  53  shows  three  of  the  possible  scenarios  that  may 
occur  if  a  second  clockwise  turn  is  added  to  the  topology  shown  in  Fig.  52(c).  Fig.  53(a)  is  illegal 
because  it  overlaps  the  initial  straight  section.  Fig.  53(b)  is  a  legal  topology  where  the  initial  straight 
section  has  been  shortened.  Fig.  53(c)  is  a  legal  topology  where  the  first  clockwise  turn  has  been 
lengthened.  As  sections  are  added,  the  problem  becomes  more  combinatorial.  We  eliminate  this 
complication  by  focusing  on  standardized  serpentine  and  spiral  topologies,  which  are  both  highly 
compact  and  readily  integrability  with  other  on-chip  components. 


The  definition  of  a  serpentine  or  spiral  channel  topology  can  be  contained  in  an  object  .  Ser¬ 
pentine  topologies  can  be  represented  as  =  \Lslr hLlrn 2,Ls,r ^,Llnl 4,.. . ,  Lstrn,Ltrnm].  The  odd  ele- 


section.  (c)  Legalize  by  lengthening  the  first  clockwise  turn. 
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FIGURE  54.  Bounding-box  determination  for  (a)  serpentine  topologies  and  (b)  spiral  topologies.  All  points 
are  calculated  with  respect  to  Pg  (the  origin).  The  X  and  Y  dimensions  are  the  actual  maximumtopology 
dimensions  plus  2PAD,  where  PAD  is  a  minimum  feature  size  specification. 

ments  of  are  straight  section  lengths  Lstrn  where  n  =  1,3,5...  |  |-1,  and  the  even  elements  are  turn 
lengths  Ltrnm  along  the  center-line  radius  of  the  channel,  where  m  =  2,4,6...  |  |.  Spiral  topologies 

can  be  represented  as  an  array  =  \Ltrn  1,Ltm2,...jJrnm ].  The  length  of  each  U-bend  Ltrnm  is  mea¬ 
sured  along  the  center  line  radius  of  the  channel  and  m  =  1 _ |  |. 

Once  a  topology  has  been  specified,  the  space  required  to  fit  it  on  a  chip  can  be  detennined 
by  measuring  the  area  or  bounding-box  of  the  topology.  Fig.  54  shows  how  the  bounding  box  of  a 
serpentine  and  spiral  topology  is  determined.  A  set  of  points  P  along  the  edges  of  the  topology  are 
detennined  starting  from  point  p0  =  (0,0),  p0  e  P.  All  subsequent  points  are  tracked  with  respect  to 

Po- 

The  width  Xand  height  7  of  a  topology  can  be  readily  determined  according  to  Eqs.  (58)  and 

(59). 

X  =  max(x/,  Vxz-  e  P)  -  min(x/,  \/xt  e  P)  (58) 

7  =  max(v,,  Vv(  eF)  -  min(v;,  Vv,  e  P)  (59) 

The  simulator  described  in  Section  4.2  as  well  as  the  connectivity  and  area  constraints  are 
used  to  evaluate  each  candidate  design  as  shown  in  Eq.  (60).  This  function,  labelled  ASIM,  is  capa- 


ble  of  simulating  and  calculating  the  bounding  box  dimensions  of  any  arbitrary  topology. 

[X,  Y,  E,  Rjj,  Cj]  =  ASIM(V,geom, props,  buffer)  (60) 

The  simulator  takes  in  the  voltage  drop  V  along  the  channel,  and  the  topology  object  geom 
where  the  instance  e  geom.  The  props  object  which  contains  analyte  physiochemical  properties 
and  the  buffer  object  which  contains  the  buffer  physiochemical  properties  are  also  passed  into  the 
simulator.  The  bounding  box  dimensions,  X  and  Y,  are  calculated  internally  and  returned  along  with 
the  electric  field  strength  E  as  well  as  the  separation  resolution  Rjj  and  peak  concentrations  Q  for 
each  band  at  the  end  of  the  channel.  ASIM  forms  the  core  of  the  NLP  fonnulations  that  we  present 
in  the  following  section  because  it  maps  the  design  variables  (channel  section  lengths  and  voltage) 
to  device  size  and  performance. 

4.3.2.  Non-Linear  Programming  Optimization 

Here  we  describe  non-linear  programming  (NLP)  optimization  approaches  and  problem  for¬ 
mulations  to  directly  optimize  the  serpentine  and  spiral  topologies  [101].  Our  approach  is  to  repeat¬ 
edly  solve  optimization  problems  for  each  topology,  incrementing  the  number  of  sections  in  the 
topology  up  to  a  value  NS,  which  is  the  maximum  number  of  channel  sections  that  can  fit  within  the 
bounding  box  defined  by  Xmax  and  Ymax.  NS  is  readily  found  using  a  channel  placement  and  pack¬ 
ing  heuristic  [102].  Optimizing  each  NLP  individually  is  effective  since  the  problem  is  not  combi¬ 
natorial. 

Our  approach  results  in  NS  locally  optimal  solutions  describing  how  topology  size  changes 
with  the  number  of  sections  in  the  topology.  There  are  several  methods  available  that  can  be  used  to 
solve  the  serpentine  and  spiral  optimization  problems  [103],  We  were  able  to  apply  a  standard  suc¬ 
cessive  quadratic  programming  (SQP)  [43]  algorithm  to  solve  the  spiral  optimization  problem. 
However,  a  more  tailored  approach  was  required  to  optimize  the  serpentine  because  our  objective 
function  is  non-smooth,  which  results  in  discontinuous  gradients.  The  non-smoothness  results  from 
the  fact  that  both  increasing  and  decreasing  the  length  of  an  inter-tum  straight  section  can  poten¬ 
tially  result  in  an  increase  in  design  area.  Furthennore,  the  problem  is  non-convex  and  the  con¬ 
straints  can  become  difficult  to  satisfy  when  certain  physical  property  values  or  bound 
specifications  are  used. 

Adaptive  Penalty  Parameters.  A  standard  penalty  function  method  with  fixed  penalties  must  be 
minimized  several  times  to  satisfy  the  constraints  of  the  original  NLP  to  a  high  level  of  accuracy.  A 


82 


common  way  to  deal  with  this  inefficiency  is  to  minimize  an  exact  penalty  function  where  the  L  / 
norm  of  the  constraint  violation  is  added  to  the  objective  function.  If  a  sufficiently  large  penalty 
parameter  is  chosen,  the  solution  to  the  original  NLP  and  the  penalized  problem  correspond  exactly; 
however,  the  penalized  problem  is  non-smooth.  We  found  that  in  our  work,  the  exact  penalty  formu¬ 
lation  would  often  fail  to  converge.  We  speculate  that  this  occurred  because  the  large  penalty 
parameters  required  for  constraint  satisfaction  also  introduce  ill  conditioning  in  the  Hessian.  Also, 
the  large  number  of  constraints  in  our  problem  introduce  a  great  deal  of  non-smoothness  into  the 
problem. 

Instead,  we  have  chosen  to  address  this  problem  by  changing  the  scalar  penalty  parameter  pk 
into  an  adaptive  penalty  parameter  pk  as  shown  in  Eq.  (61)  for  each  of  the  constraints  /  =  1  ...  in  in 

the  original  NLP.  pk  is  calculated  based  on  the  scaled  constraint  violation  from  the  previous  itera¬ 
tion  k- 1 . 


min  F(xk,  pk)  =  X  ■  Y  +  V  pf 


Lx  A 


i  =  l ... m 


^max(0,g(/  ' 


max(0,  g,(x  )) 
typX 


(61) 


typX 


Vi 


1 . . .  m 


InEq.  (61),  xk  =  {  ,V}  for  the  current  iteration,  p  is  an  increasing  sequence  10,100,1000,  ..., 
and  typX  is  a  typical  x  value  meant  to  scale  g,(x)  so  that  all  constraint  residuals  are  of  roughly  the 
same  order. 

This  is  an  intuitive  approach  based  on  the  fact  that  the  Lagrange  multiplier  Xk  «  pk~ 
max(0,g(x^))  and  it  approaches  its  optimal  value  X*  as  F  and  x  approach  their  optimal  values  F*  and 


x*.  The  scalar  parameter  pk  is  necessary  to  ensure  that pk  >  Xk  so  that  each  constraint  is  appropri¬ 
ately  enforced. 

The  procedure  begins  by  obtaining  p/  from  the  initial  point  x°.  Subsequent  values  of pk  are 


calculated  based  on  the  solution  to  the  previous  iteration  xk~] .  Our  computational  experience  indi¬ 
cates  that  our  adaptive  penalty  parameter  approach  typically  requires  1  to  2  fewer  sequential  mini¬ 
mizations  than  the  standard  approach. 
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Table  13. Spiral  Design  Data  (extracted  from  [20]). 


Data _ Values 

Actual  Design  Area 


Total  Separation  Length 
Channel  Width 
Channel  Depth 
Applied  Electric  Field 
Applied  Voltage  (over  22.2  cm) 
Maximum  Available  Voltage 
NUmber  of  Plates  (@22.2  cm) 


3.7  cm  x  3.8  cm  (A  =  14  cur) 

22.2  cm 

40  pm 

15  pm 

1170  Vein1 

26.0  kV 

30.0  kV 

NDCF  =  1.04  x  106 _ 


4.3.3.  Example  Case  Studies 

Here  we  compare  designs  generated  using  our  techniques  to  a  serpentine  design  developed 
by  Jacobson  et  al.  [62]  and  a  spiral  design  developed  by  Culbertson  et  al.  [20].  We  attempted  to 
include  industrially  relevant  fabrication  specifications.  An  minimum  feature  size  specification  of 
PAD  >  5  pm  was  obtained  from  a  MEMS  and  micro  fluidics  foundry  [104].  Using  our  methods,  we 
were  able  to  reduce  the  size  of  the  literature  designs  while  still  maintaining  the  desired  level  of  sep¬ 
aration  performance. 

Optimization  of  a  Spiral  Design.  Culbertson  et  al.  present  a  spiral  design  [20]  with  4  channel  sec¬ 
tions  between  the  injector  and  detector.  The  entire  design,  including  all  wells,  fits  within  a  5cm  x 
5cm  area.  The  spiral  topology  itself  fits  within  approximately  a  3.7cm  x  3.8cm  area.  In  our  compar¬ 
ison  we  will  focus  only  on  the  area  occupied  by  the  spiral  topology.  The  chemical  species  used  was 
dichlorofluoroscein  (DCF)  and  its  physical  properties  were  extracted  from  results  presented  by  Cul¬ 
bertson  et  al.  The  buffer  used  was  a  boric  acid/TRIS  solution.  Table  13  lists  some  of  the  important 
design  features  and  data. 

Fig.  55(a)  shows  a  schematic  of  the  separation  portion  of  the  original  spiral  topology.  Here 
we  use  the  optimization  methods  described  in  section  Section  4.3.2  to  design  a  spiral  separation 
system  that  is  smaller  than  the  original  design  and  maintains  the  same  or  greater  plate  number 

(Ndcf>  1.04  x  106). 

Fig.  55(a),  shows  the  original  design,  which  is  essentially  4  semicircles  with  large  radii  and 
narrow  channel  widths.  This  configuration  leads  to  low  dispersion.  However,  the  design  leaves  a 
great  deal  of  inter-channel  spacing  which  allows  us  to  significantly  compact  the  design.  In 
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X(pm)  x104  X(|im)  x1Q4 

(a)  (b) 

FIGURE  55.  Improvement  of  a  spiral  design:  (a)  A  schematic  of  the  design  presented  in  |20]  (b)  A  schematic 
of  an  improved  design  generated  using  our  optimization  algorithm. 

Fig.  55(b)  we  see  that  the  design  can  actually  be  reduced  in  area  by  approximately  56  %. 

This  results  from  the  addition  of  6  new  channel  sections  resulting  in  a  10  section  design.  An 
applied  voltage  of  27  kV  was  used,  which  leaves  an  extra  3  kV  to  connect  the  separation  channel  to 
the  wells. 

The  channels  are  compressed  as  closely  as  possible  to  each  other  without  violating  the  con¬ 
straint  on  PAD  >  5  pm.  Fig.  55  shows  a  magnification  of  the  optimized  design.  It  can  be  seen  that 

the  feature  spacing  requirement  is  met.  Furthermore,  this  design  retains  approximately  4.5  cm  of 


FIGURE  56.  Magnification  of  the  optimized  spiral  design  shown  in  Fig.  55. 
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space  interior  to  the  spiral  for  the  placement  of  the  waste  well.  This  design  could  be  further  com¬ 
pressed  by  either  adding  more  sections,  increasing  the  applied  voltage  or  decreasing  the  required 
inter-channel  spacing. 

4.3.4.  Summary 

In  this  section  we  formulated  the  layout  optimization  problem  for  serpentine  and  spiral 
topologies.  We  investigated  numerical  techniques  to  directly  solve  our  NLP  formulations.  We  found 
that  while  a  conventional  SQP  algorithm  is  effective  for  spiral  topologies,  other  methods  must  be 
applied  for  serpentine  topologies. 

4.4.  Place  and  Route 

Multiplexed  LoCs  have  many  potential  areas  of  application,  especially  in  the  biomedical 
industry,  where  device  size  is  highly  constrained.  For  instance,  in-vivo  diagnostic  chips  that  are 
implanted  into  living  tissue  must  be  both  small  and  contain  a  high  degree  of  functionality.  Hand¬ 
held  point-of-care  (PoC)  medical  devices  are  another  important  example.  In  these  devices,  size  is 
not  only  a  major  constraint,  but  also  the  primary  cost  consideration,  since  LoCs  for  PoC  applica¬ 
tions  are  often  disposable. 

Currently,  the  most  complex  multiplexed  LoCs  consist  of  arrays,  where  hundreds  of  repli¬ 
cated  simple  channel  structures  function  in  parallel  [3].  The  typical  layout  of  these  devices  is  in  a 
spoked-wheel  configuration  with  straight  separation  channel  sections  making  up  the  spokes.  Each 
subsystem  design  is  defined  by  the  set  of  chemical  properties  and  design  specifications  discussed  in 
Section  4.3. 

4.4.1.  Floorplan  Problem  Definition 

We  define  our  floor  plan  problem  as  follows  [105]:  given  a  set  of  subsystems  N  and  their 
associated  input  and  output  wells  obtain  a  planar  arrangement  A  of  subsystems  and  wells  that 
does  not  exceed  the  total  chip  height  H  or  width  W,  and  where  all  subsystems  with  respect  to  the 
constraint  set  in  the  serpentine  optimization  problem  are  feasible.  The  objective  is  to  choose  the 
arrangement  of  subsystems  and  wells  that  is  most  compact  and  where  fluid  I/O  wells  are  positioned 
on  the  edge  closest  to  their  associated  subsystem  ports. 

We  divide  the  fonnulation  into  5  parts:  PI  -  Subsystem  rotation,  P2  -  Subsystem  overlap,  P3 
-  Well  placement,  P4  -  Well  to  port  assignment,  and  P5  -  Boundary  constraints  and  objective. 
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FIGURE  57.  Eight  possible  subsystem  orientations. 


Subsystem  Orientation  (PI).  In  Eq.  (62)  we  define  the  8  possible  orientations  and  associated  port 
locations  for  each  subsystem  i  in  K  The  8  orientations  are  rotations  and  reflections  of  a  serpentine 
subsystem  as  shown  in  Fig.  57. 
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Vr  =  1 ...  8 


4 

St  =  (S\  S?) 
Sw{  =  ( Swx ,  Swy) 
B,  =  (B\ By) 
Bw-  =  ( Bwx ,  By/) 


Vi  €  K 


(62) 


In  Eq.  (62),  one  of  the  Boolean  variables,  Z  ;  through  Z  ,  will  be  assigned  True  depending  on 
it's  orientation  as  shown  in  Fig.  57.  When  Boolean  variables  are  translated  into  mixed  integer  form, 
they  are  assigned  a  binary  value,  True  =  1  or  False  =  0. 

The  coordinates  of  the  sample  (5)),  sample  waste  (Sw,),  buffer  (Bx),  and  buffer  waste  (Bw,) 
ports  are  determined  for  an  orientation  of  a  subsystem  with  respect  to  the  bottom  left-hand  comer  of 
a  particular  block.  For  example,  the  buffer  port  location  Bx  when  Z4 { =  True  is  Bi  =  (Bx,  Bv)  =  (xx  +  Yx 
-co  -  PAD,  Vj  +  Xj).  The  port  locations  are  a  function  of  the  subsystem’s  location,  x,  and  yx,  length  X-t 
and  height  Yn  as  well  as  the  minimum  feature  spacing  PAD  and  channel  width  co. 


Overlap  Prevention  (P2).  Since  we  assume  a  building  block  layout  style,  blocks  can  be  placed 
anywhere  on  the  chip,  but  a  legal  placement  is  only  achieved  when  no  blocks  overlap.  Overlap  is 
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prevented  by  assigning  {right,  left,  above,  below }  relationships  between  each  block  in  K , 
(i,j  e  K  ).  Here  we  show  a  GDP  that  enforces  these  relationships  for  i  left  of  j  (Eq.  (63)),  /'  right  of 
j  (Eq.  (64)),  i  below  j  (Eq.  (65))  and  i  above  j  (Eq.  (66)). 
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ord(i)  <  ord(j),  \/i,j  e  K 

The  values  of  the  binary  variables  81  tj,  ir  LJ,  83 \j,  and  li4 LJ  detennine  the  orthogonal  packing 
of  the  blocks.  A  disjunction  is  active  when  its  associated  binary  equals  one.  The  disjunctions  over 
Zf  incorporate  the  block  orientation  information  into  the  block  packing  problem. 


Well  Placement  (P3).  We  require  that  fluid  I/O  wells  are  placed  on  the  chip's  periphery  as  shown  in 
Fig.  58.  Eq.  (67)  shows  the  nested  disjunctions  necessary  to  model  the  placement  of  a  well  on  the 
top  edge  of  the  chip.  We  use  the  variables  cph,  cpE,  g>hi  and  to  determine  if  well  ie  is  placed  on 
the  top,  right,  bottom  or  left  edge  of  the  chip  respectively.  The  value  of  (j)y  dictates  the  relative  posi- 


88 


FIGURE  58.  (a)  Schematic  of  a  serpentine  subsystem  showing  horizontal  JV  and  vertical  V  dimensions,  and 
port  locations  (.S'A,.S'-'  ),  (Bx,By),  (Swx^wy),  and  ( B\P,Bw y).  (b)  Three  subsystem  chip  layout  showing  the 
horizontal  W  and  vertical  H  dimensions  and  subsystem  locations  (x^y);.  Fluid  I/O  wells  are  shown  as 
circles,  and  auxiliary  channels  (black  lines)  connect  wells  and  subsystems. 

tion  of  well  i  with  respect  to  well  j  e  if  it  exists  on  that  particular  edge.  A  well's  position  is 

defined  by  a  coordinate  point  in  the  center  of  the  well  Its  placement  is  a  function  of  the 

the  chip's  width  W  and  height  H,  as  well  as  the  well's  diameter  d  and  the  PAD  around  every  well. 

The  disjunctions  for  the  right,  bottom,  and  left  edges  are  similar  to  Eq.  (67). 
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linearity  to  the  system. 

Well  to  Port  Assignment  (P4).  Here  we  model  the  connection  between  subsystem  ports  and  wells 
as  an  assignment  problem.  We  require  that  there  is  a  unique  one  to  one  mapping  between  ports  and 
wells.  Our  goal  is  to  minimize  the  distance  that  auxiliary  channels  must  span  to  connect  a  port  to  a 
well.  Eqs.  (68)  enforce  that  a  port  of  a  particular  subsystem  is  assigned  to  only  one  well. 


I  Pj  = 1 

JeW 


Z  Pj  =  1 

jeW 


V7  €  N 


Z  = 1  Z  PJ  - 1 

j  e  W  j  eW 


(68) 


The  binary  variables,  PSjj,  PSwij,  PBjj  and  PBwjj  represent  the  connection  of  the  sample, 
sample  waste,  buffer  and  buffer  waste  ports  to  a  particular  well.  Eqs.  (69)  enforces  that  a  well  is 
connected  to  a  single  port. 


,,S’w 


y  (p° . + pow  +  pB .  +  pBw) 
zL  v  lJ  l’J  J 


V/  e  W 


(69) 


We  use  a  rectilinear  distance  metric  [106]  to  estimate  the  length  that  an  auxiliary  channel 
must  span  between  a  port  z  e  N  and  a  well  /  e  W .  The  connection  length  between  a  particular  port 
and  well  represents  the  cost  of  an  assignment  and  is  captured  by  the  variables  CSj  p  CSwjj,  CB LJ  and 

CBwjj  as  shown  in  Eq.  (70).  Notice  that  the  cost  of  an  assignment  is  not  fixed  as  it  is  in  typical 
assignment  problems,  but  is  a  function  of  the  chip's  size  and  the  size  and  position  of  each  sub¬ 
system. 
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(70) 


In  practice,  Eqs.  (70)  are  implemented  as  a  system  of  linear  inequalities  which  avoids  the 
nonlinearity  introduced  by  the  abs{ •)  operator. 


Boundary  Constraints  and  Objective  (P5).  The  bounding  box  of  the  design  is  calculated  in  Eqs. 
(71),  where  W  and  H  are  the  calculated  design  width  and  height  respectively.  Notice  that  all  sub- 
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systems  and  wells  are  taken  into  account. 


4  4 

xt  +  Xt  £  Z{r~  1  +  7.  ^  Z(2r  +  cp"(J  +  2PAD)  +  (pj(f/  +  2PAD)  <  W 

(71) 

4  4 

y(.  +  7.  ^  Z,2r"  1  +  X.  ^  Z(2r  +  cp5(rf  +  2PAD)  +  cpf  ( d  +  2PAD)  <  W 

r—  1  r =  1 

When  the  available  chip  width  IF  and  height  PI  are  specified,  IF  and  PI  must  be  constrained 

such  that  W  <W  and  H  <H .  However,  these  constraints  add  a  great  deal  of  computational  burden 
to  an  already  difficult  problem. 

Here  we  have  chosen  to  use  the  weighted  sum  of  the  bounding  box  perimeter  and  the  total 
auxiliary  channel  length  as  our  objective  function  (Eq.  (72)).  The  weights,  a,;  and  a2,  are  assigned 
values  to  emphasize  either  the  area  minimization  or  the  auxiliary  channel  length  minimization. 
However,  choosing  good  values  for  a  1  and  a2  is  often  difficult.  We  will  discuss  our  solution  to  this 
issue  later  in  the  report. 


mini7  =  al(W+H)  +  a2  j  +  CUJPUJ  +  cl/lj  +  cf,  P^J  )  (72) 

i  <E  t\  j  <E  W 


Problems  PI  -  P5  represent  a  rigorous  yet  highly  coupled,  highly  combinatorial  nonlinear 
description  of  the  multiplexed  LoC  design  and  layout  problem.  This  kind  of  fonnulation  is  typically 
translated  into  a  mixed  integer  nonlinear  program  (MINLP)  and  solved  using  detenninistic  optimi¬ 
zation  approaches  that  combine  gradient-based  search  methods  with  Branch  and  Bound  [107]. 
However,  the  floor  planning  problem  that  we  have  formulated  is  a  rectangle  packing  problem  that 
has  been  shown  to  be  NP-hard  [108],  In  addition,  floor  plan  optimization  using  Branch  and  Bound 
[109]  or  conventional  mixed  integer  linear  programming  (MILP)  approaches  [110]  often  have  diffi¬ 
culty  handling  problems  of  realistic  size.  Since  we  are  coupling  rectangle  packing  with  a  non-con- 
vex  (nonlinear)  physiochemical  problem,  our  problem  is  more  difficult  than  rectangle  packing 
alone.  In  fact,  even  if  we  neglect  parts  PI,  P3  and  P4  of  the  problem,  we  can  only  solve  for  trivial 
instances  (i.e.  X  <  4  and  W  =  0)  using  conventional  MINLP  solvers. 


4.4.2.  Combinatorial  Floorplan  Problem  Formulation 

We  now  formulate  the  multiplexed  LoC  floorplanning  problem  as  a  combinatorial  optimiza- 
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(a)  (b)  (c) 

FIGURE  59.  (a)  One  serpentine  unit  (x  =  1).  (b)  Two  serpentine  units  (x  =  2).  (c)  Completing  the  topology 
with  an  initial  straight  section. 

tion  problem  directly  (the  above  GDP  approach  clearly  requires  a  heuristic  for  solution). 

Subsystem  Orientation.  As  discussed  earlier,  each  subsystem  can  have  one  of  eight  possible  orien¬ 
tations  (Fig.  57).  Different  orientations  contribute  to  design  compaction.  We  define  a  vector 

Z  =  (zl5z2, _ , Z|K|)  which  contains  the  subsystem  orientation  labels  z,-  e  {1,  ...,  8}  where  i  is  a 

unique  element  of  K  .  Each  z;  defines  the  orientation  of  subsystem  i  and  allows  for  the  deterministic 
calculation  of  that  subsystem's  port  locations.  For  example,  the  buffer  port  location  Bl  for  orienta¬ 
tion  Zj  =  4  is  Bj  =  (Bx,  Bv)  =  (Xj  +  Wj  -  co  -  PAD,  yi  +  hj),  which  corresponds  to  the  case  when  z,4  = 
True  in  the  GDP.  The  port  locations  are  a  function  of  the  subsystem's  location  x;  and  y,-,  width  w-t  and 
height  hi,  as  well  as  the  feature  spacing  PAD  and  channel  width  co.  However,  recall  that  the  perfor¬ 
mance  simulation  has  no  knowledge  of  subsystem  orientation.  Thus,  Eq.  (73)  must  be  applied  to 
translate  the  general  subsystem  dimensions,  Xt  and  Yh  used  in  the  simulator  to  the  actual  subsystem 
dimensions  vv;  and  hr 

\  xi  if  odd 

wi  =  ) 

[  T-  otherwise. 

Serpentine  Topology  Size.  The  types  list  can  be  reduced  to  a  scalar  value  x  indicating  the  number 
of  serpentine  units  within  a  subsystem  topology  as  shown  in  Fig.  59.  Fig.  59(a)  and  (b)  show  a 
topology  with  x  =  1  and  x  =  2  respectively.  The  topology  is  completed  by  adding  an  initial  straight 


f  Xi  if  z  ■  even 

hi  =  1  (73) 

( 7-  otherwise. 


$ 

section  to  the  serpentine  as  shown  in  Fig.  59(c).  We  define  a  vector  S  =  (x1?  x2,  ...,  X|K|)  which 

UB 

contains  the  number  of  serpentine  units  x;  e  {  1,  ...,  x  }  for  each  subsystem,  i  e  K  .  Here,  x,-  can 

range  from  1  to  a  user  specified  upper  bound,  xUB.  Our  experience  indicates  that  a  conservative 

value  for  xUB  is  15,  which  corresponds  to  a  serpentine  composed  of  61  individual  channel  sections 
where  the  number  of  sections  =  4-  x(-  +  1 . 


Relative  Subsystem  Positions.  We  encode  the  constraints  to  prevent  subsystems  from  overlapping 
using  a  VLSI  floor  plan  representation  known  as  the  Sequence  Pair  (SP)  [108].  SP  encodes  the  {left, 
right,  above,  below }  spatial  relations  between  objects  in  the  plane  into  a  concise  structure.  This 
structure  is  readily  searchable  using  standard  heuristic  methods  such  as  Simulated  Annealing  (SA) 

[HI]- 

A  SP  is  composed  of  two  K  -tuples  T  =  (r+,r_),  where  T+  and  T_  are  ordered  lists  of  the  sub¬ 
system  labels  in  K  .  A  SP  maps  to  |K|(|K[—  l)/2  linear  math  programming  constraints  as  shown 
in  Eq.  (74). 


xi  +  wi<xj<^(...i...j...,  ...) 

yi  +  hi  .../... y...) 


(74) 


Here,  i  and  j  are  unique  subsystem  labels  in  K  .  If  i  appears  before  j  in  both  T+  and  T_,  then 
subsystem  i  is  left  of  j.  Likewise,  if  i  appears  after  j  in  T+  and  before  j  in  T_,  then  subsystem  i  is 
below  j.  The  SP  representation  has  several  attractive  features.  First,  the  solution  space,  although 

large,  is  finite  i.e.  cc  (|K|!)  .  Second,  the  neighborhood  of  a  particular  arrangement  A  is  readily 
constructable  through  simple  perturbations  of  L  Furthermore,  every  evaluation  of  T  results  in  a  fea¬ 
sible  planar  arrangement,  which  is  not  true  for  all  non-slicing  floorplan  representations  (e.g.  Comer 
Block  List  presented  by  Hong  et  al.  [112]).  Finally,  SP  is  a  general  floorplan  representation,  which 
means  that  optimal  solutions  are  not  excluded  due  to  assumptions  about  floorplan  structure. 


Well  Placement.  Confining  the  fluid  I/O  wells  to  the  edges  of  the  chip  facilitates  fluid  transport  on 
and  off  the  chip.  We  have  developed  a  graph-based  approach  that  places  a  well  on  the  chip  edge 
closest  to  its  associated  port. 

For  a  given  arrangement  A  of  subsystems,  we  construct  an  undirected,  edge- weighted  planar 
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(a)  (b) 

FIGURE  60.  (a)  Compacted  arrangement  of  subsystems,  (b)  Corresponding  well  placement  graph. 

graph  G,j  =  <V,E>  by  noting  that  a  compacted  arrangement  of  subsystems  is  essentially  a  planar 
graph  (Fig.  60).  The  graph  we  create  is  similar  to  a  placement  graph  used  in  VLSI  circuit  design 

[113]. 

We  construct  the  graph  by  first  generating  a  set  of  preliminary  nodes  by  taking  the  Cartesian 
product  of  the  bounding  box  dimensions,  subsystem  locations  and  port  locations  as  shown  in  Eq. 
(75). 


A  set  of  preliminary  edges  are  fonned  by  connecting  each  node  with  its  nearest  neighbor  in 
the  x  and  y  directions.  The  graph's  vertex  set  V  and  edge  set  E  are  generated  by  removing  all  the 

nodes  and  associated  edges  that  reside  within  subsystem  boundaries.  Finally,  four  supernodes,  v^, 

vT,  vR,  and  vB  are  added  to  the  vertex  set.  The  supernodes  are  connected  to  the  left-most,  top-most, 
right-most,  and  bottom-most  nodes  of  GA  as  shown  in  Fig.  60(b). 

The  weights  assigned  to  edges  within  the  interior  of  GA  are  simply  the  distances  in  the  jc  and 
y  directions  between  connected  adjacent  nodes.  The  edges  connecting  supernodes  are  weighted 
based  on  the  maximum  bounding  box  dimension.  This  prevents  short-circuit  paths  through  the 
graph  during  the  well  placement  procedure.  Each  well  is  placed  on  one  of  the  four  sides  of  A  as  fol¬ 
lows: 
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(a)  (b) 

FIGURE  61.  (a)  Poor  well  placement  based  on  rectilinear  model  (dotted  arrow).  (b)  Compacted 
arrangement  pulled  apart  by  weighted-sum  objective  function. 


1.  Construct  G/t  and  initialize  sets:  left  =  right  =  top  =  bottom  =  0. 

2.  Select  a  port  vertex  p  e  vp,  where  vpa  V  is  the  set  of  ports. 

3.  Determine  the  shortest  paths  from  p  to  each  supemode  vL,  vT,  vR,  v8.  The  shortest 
distance  among  these  detennines  a  connection  from  p  to  that  edge.  Store  this  distance 
in  the  set  P. 

4.  Assign  p  to  left  if  p  connects  to  vL,  to  top  if  p  connects  to  vT,  to  right  if  p  connects  to 
vR  or  to  bottom  if  p  connects  to  vB. 

5.  Order  wells  in  left  or  right  based  on  the  v-coordinates  of  their  ports. 

6.  Order  wells  in  top  or  bottom  based  on  the  ^-coordinates  of  their  ports. 

7.  Obtain  a  total  routing  length  estimate,  (j),  by  summing  the  entries  in  P.  (note  that 
\P\  =  \W\  =  4-  |K|) 

In  step  3,  we  apply  a  shortest  path  algorithm  to  efficiently  search  GA.  Steps  5  and  6  can  be 
performed  using  an  efficient  sorting  method. 

Our  solution  methodology  allows  us  to  avoid  two  fundamental  problems  that  would  result 
from  a  more  simplistic  formulation.  The  first  possible  problem  is  illustrated  in  Fig.  61(a).  Here 
wells  are  placed  using  a  standard  rectilinear  distance  metric.  However,  the  shortest  rectilinear  dis¬ 
tance  between  a  port  and  a  well,  indicated  by  the  dotted  arrow,  may  pass  through  a  subsystem  and  is 
therefore  a  poor  estimate  of  the  true  distance  between  a  port  and  the  chip  edge. 

The  distance  estimates  produced  using  GA  are  far  more  accurate  because  they  account  for 
the  fact  that  routes  must  go  around  subsystems. 

A  second  problem  can  arise  if  a  common  VLSI  floorplanning  objective  function  of  the 
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Boundary  constraints 
Subsystem  performance 
Subsystem  operation 
Fabrication  constraints 


FIGURE  62.  Flowchart  for  Floorplan  Solution  Method. 


form:  a,  •  Area  +  a2  •  WL  is  used,  where  Area  is  chip  area,  WL  is  wire  length  and  o.|  and  a2  are 

weighting  factors.  This  objective  works  well  for  typical  circuit  design  problems  because  circuit  ele¬ 
ments  are  highly  interconnected  (i.e.  wired  together)  within  the  chip  interior  and  a  wide  variety  of 
a  |  and  a2  values  will  produce  satisfactory  results.  Since  we  connect  single  ports  within  the  chip's 
interior  to  single  wells  on  the  chip's  edges,  over-weighting  a2  pulls-apart  a  compact  arrangement 
and  results  in  a  large  amount  of  unused  space  within  the  chip  interior  (Fig.  61(b)).  Since  we  gener¬ 
ally  do  not  know  good  values  for  a.\  and  a2  a  priori,  we  first  compact  the  subsystems  and  lock  their 
positions.  Then  we  place  the  wells,  thereby  decoupling  the  continuous  variables  of  the  problem 

while  retaining  the  global  nature  of  the  combinatoric  variables,  T,  Z  and  S . 

Floorplan  Solution.  A  flowchart  for  the  floorplanning  algorithm  is  shown  in  Fig.  62.  The  main 
idea  is  to  use  a  probabilistic  search  heuristic  such  as  SA  to  deal  with  the  combinatorial  aspects  of 
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the  problem  and  an  efficient  gradient  based  method  for  the  remaining  continuous  space  problem. 
We  have  chosen  this  hybrid  approach  because  SA  has  been  shown  to  be  an  effective  search  method 
for  difficult  combinatorial  problems,  but  is  generally  inferior  to  gradient-based  methods  on  well 
posed  continuous-space  problems  [114]. 

4.4.3.  Routing  of  Auxiliary  Microfluidic  Channels 

The  goal  of  the  routing  stage  is  to  take  a  compacted  arrangement  generated  by  our  floorplan¬ 
ning  algorithm  and  to  determine  the  final  placement  or  precise  positions  of  subsystems  and  auxil¬ 
iary  channels  that  connect  ports  to  wells.  The  auxiliary  channels  are  kept  as  short  and  straight  as 
possible  to  reduce  fabrication  costs  and  to  allow  for  more  convenient  fluid  loading  and  dispensing. 
In  addition,  short  straight  channels  reduce  dispersion  effects  [115]  and  concentration  non-unifor¬ 
mity  [26]  as  well  as  the  operating  voltage  required  to  drive  the  chip. 

Routing  Problem  Description.  We  define  our  routing  problem  as  follows:  Given  a  compact 
arrangement  A  of  subsystems,  find  a  routing  solution  which  consists  of  the  set  of  planar  paths  P 
through  A  such  that  each  port  is  connected  to  a  single  well  and  that  the  total  length  and  number  of 
bends  in  P  is  minimized. 

The  routing  of  auxiliary  channels  is  similar  to  wire  routing  in  VLSI  circuit  design,  however, 
there  are  several  complicating  features.  Most  importantly,  LoC  devices  are  generally  fabricated  in  a 
single  layer  to  reduce  fabrication  costs  and  complexity  [115].  This  means  that  all  auxiliary  channels 
must  be  routed  in  a  planar  fashion  and  can  not  be  routed  above  or  below  a  subsystem.  Furthermore, 
the  assumptions  that  channels  can  feed  through  a  subsystem  or  that  ports  may  move  along  a  sub¬ 
system  edge  [113],  do  not  apply  in  our  fonnulation.  Additionally,  auxiliary  channels  occupy  signif¬ 
icant  space  on  the  chip,  and  therefore,  the  exact  placement  of  routes  is  critical  to  the  quality  of  the 
overall  design. 

The  single-layer  routing  problem  has  been  investigated  in  the  VLSI  literature  [113].  How¬ 
ever,  much  of  this  work  is  not  directly  relevant  to  our  problem  because  it  either  pertains  to  global, 
multi-terminal,  or  detailed  routing  in  bounded  regions  (e.g..  channel  routing).  In  global  routing,  the 
paths  that  wires  must  follow  around  obstacles  on  the  chip  are  detennined  in  a  general  way.  A  subse¬ 
quent  detailed  routing  stage  places  the  wires  precisely.  Also,  most  VLSI  routing  algorithms  con¬ 
struct  multi-terminal  routes,  since  single  wires  must  often  connect  several  terminals  on  the  chip. 
Unlike  channel-routing,  which  takes  place  in  small  bounded  regions  of  the  chip,  we  are  interested  in 
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finding  routing  solutions  through  the  whole  chip.  Finally,  many  general  VLSI  routing  algorithms 
employ  a  sequential  rip-up  and  re-route  approach  to  construct  a  routing  solution.  This  approach 
lacks  feasibility  guarantees  and  can  become  artificially  over-constrained  when  poor  routing  choices 
are  made  early  in  the  procedure. 

We  are  currently  developing  a  routing  procedure  that  is  specifically  designed  to  solve  our 
two-terminal,  single-layer,  detailed  routing  problem.  In  addition,  we  directly  address  length  and 
bend  minimization.  In  our  approach,  we  hope  to  exploit  the  network  integrality  property  of  the  min¬ 
imum-cost  network  flow  (MCNF)  fonnulation  to  find  routing  solutions  in  a  simultaneous  fashion 
using  linear  programming  (LP).  Large-scale  MCNF  problems  can  be  solved  using  standard  LP  solv¬ 
ers  or  tailored  algorithms  [116].  Furthennore,  LP  can  quickly  provide  rigorous  feasibility  guaran¬ 
tees.  However,  the  bend-minimization  constraints  discussed  in  the  following  section  break  the 
MCNF  structure,  thus  forcing  us  to  use  an  integer  programming  (IP)  approach.  While  many  differ¬ 
ent  IP  formulations  are  possible  to  solve  our  routing  problem,  it  is  unclear  which  formulation  will 
perform  the  best  in  all  cases  since  there  are  non-trivial  trade-offs  between  the  number  of  constraints, 
the  type  of  constraints  and  the  number  of  variables. 

We  discuss  our  formulation,  which  is  based  loosely  on  the  concept  of  MCNF,  in  the  follow¬ 
ing  section.  We  have  successfully  used  this  formulation  to  obtain  good  solutions  for  several  small, 
compact  test  cases  and  are  currently  exploring  algorithmic  improvements  that  will  allow  us  to  deal 
with  larger  test  cases. 

4.4.4.  Routing  Solution  Method 

Fig.  63  shows  a  flowchart  of  our  procedure  for  routing  compacted  arrangements.  First,  an 
initial  arrangement  is  legalized  so  that  it  can  be  embedded  in  Ggrjcj.  Next  we  construct  Gnet  from 
Ggrid  and  use  it  to  generate  the  associated  integer  program  (IP).  The  IP  is  solved  using  a  standard 
commercial  solver  [117].  If  the  IP  is  feasible,  then  the  IP  solver  searches  for  an  optimal  routing.  If 
an  optimal  routing  or  a  routing  that  meets  a  user  specified  tolerance  is  found,  then  the  algorithm  ter¬ 
minates  successfully.  If  the  IP  is  infeasible,  the  design  is  expanded  to  open  up  new  routing  paths  and 
is  searched  again.  If  the  expanded  design  exceeds  the  available  chip  area,  the  floor  plan  algorithm 
must  be  queried  for  a  new  arrangement  that  may  have  better  routing  characteristics.  Two  of  the  key 
steps  of  our  routing  procedure  are  discussed  in  the  following  sections. 
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FIGURE  63.  Flowchart  for  the  routing  algorithm. 

Floorplan  Legalization.  Before  an  arrangement  can  be  placed  or  embedded  into  Ggricj,  it  must 
properly  confonn  to  the  minimum  allowable  feature  spacing  requirement.  If  a  subsystem's  dimen¬ 
sions  are  not  integer  multiples  of  PAD  +  co,  it  cannot  be  properly  embedded  into  the  routing  grid. 
Our  floorplan  algorithm  does  not  guarantee  that  an  arrangement  will  be  embeddable.  To  account  for 
this,  we  expand  each  subsystem  in  the  x  and  y  directions  to  the  nearest  integer  multiple  of  PAD  +  co. 
For  example,  wt  is  updated  as  follows:  w-  =  ^w/(PAD  +  Q)~|  •  (PAD  +  Q) .  The  subsystem 
height  hi  is  treated  in  a  similar  fashion.  The  current  T  value  of  the  floorplan  is  used  to  maintain  the 
relative  position  of  each  subsystem.  This  procedure  has  a  negligible  effect  on  the  perfonnance  and 
operation  of  each  subsystem  since  the  size-perturbation  represents  a  relatively  tiny  fraction  of  the 
overall  subsystem’s  size  (typically  <  1 .0%). 

The  wells  on  the  chip's  edges  are  updated  in  a  similar  fashion  so  that  they  too  conform  to 
valid  grid  points.  As  with  other  features  on  the  chip,  well  edges  are  constrained  to  be  no  less  than 
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PAD  apart.  After  legalization,  port,  well  and  subsystem  locations  confonn  to  verticies  in  Ggrj(j. 
Therefore,  when  a  routing  solution  is  found,  it  represents  the  exact  location  of  auxiliary  channels  on 
the  chip. 

Floorplan  Expansion.  For  a  given  Gnet,  the  solution  to  the  IP  formulation  will  either  yield  a  feasi¬ 
ble  routing  solution,  or  it  will  declare  the  current  routing  problem  infeasible.  If  the  IP  is  infeasible, 
then  we  are  guaranteed  that  no  routing  solution  exists  for  the  current  instance  of  Gnet. 

We  have  developed  a  procedure  that  iteratively  augments  Ggrid  by  expanding  the  arrange¬ 
ment  to  increase  the  number  of  possible  routing  paths  in  the  grid.  Each  expansion  adds  one  addi¬ 
tional  routing  path  between  each  subsystem. 

An  arrangement  is  expanded  by  first  conceptually  expanding  each  subsystem’s  width  and 
height  by  y (PAD  +  co).  We  define  the  arrangement  expansion  number  y  as  the  number  of  times  the 
arrangement  is  expanded.  The  arrangement's  T  is  used  to  maintain  the  relative  positions  of  each 
subsystem.  Finally,  each  subsystem  is  shifted  right  and  up  from  the  origin  by  y (PAD  +  co).  After 
each  expansion,  the  arrangement  is  re-placed  or  re-embedded  and  a  new  Gnet  is  generated.  Expan¬ 
sions  continue  until  either  a  feasible  routing  solution  is  obtained  or  until  the  placed  design  has 
become  too  large.  Typically,  arrangements  require  fewer  than  5  expansions  to  achieve  a  feasible 
routing  solution. 

4.4.5.  Multiplex  Synthesis  Results 

Here  we  discuss  some  of  the  key  features  and  algorithmic  behavior  of  the  floorplanning  and 
routing  algorithms.  The  results  are  encouraging  with  respect  to  both  computation  time  and  solution 
quality,  since  our  implementations  have  not  yet  been  optimized.  However,  these  results  also  indi¬ 
cate  important  areas  for  improvement  and  future  study. 

The  result  of  our  floorplanning  algorithm  for  our  5  subsystem  example  is  shown  in 
Fig.  64(a).  This  result  was  obtained  in  approximately  5  minutes  of  CPU  time  [118].  At  this  point, 
the  design  is  compact  and  all  constraints  are  completely  satisfied. 

The  result  of  the  routing  algorithm  for  our  5  subsystem  example  is  shown  in  Fig.  64(b).  This 
design  was  routed  in  approximately  3  minutes  of  CPU  time  [1 18].  At  this  point,  our  design  could  be 
fabricated. 

Floorplan  Results.  Fig.  65(a)  shows  a  placed  and  compacted  design  including  well  positions  for  a 
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FIGURE  64.  Five-subsystem  example:  (a)  Compacted  arrangement  (b)  Final  design  (routed). 

multiplexed  chip  containing  10  subsystems.  The  completed  design  shown  in  Fig.  65(b),  was  gener¬ 
ated  in  2  hours  on  a  standard  PC  [118]  with  equal  time  allocated  to  the  floorplanning  and  routing 
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Number  of  Subsystems 

FIGURE  66.  Scaling  of  NLP  compaction  subproblem  in  the  floorplanning  algorithm  per  size  of  design 
instance.  Tests  were  performed  on  a  standard  PC  (Intel  Pentium  3  CPU:  1GHz  CPU,  1Gb  RAM). 

stages.  While  we  are  pleased  with  these  results,  significantly  more  time  was  required  to  generate  the 
10  subsystem  design  than  was  required  for  the  5  subsystem  design.  The  compaction  NLP  presents 
the  most  significant  computational  bottleneck  during  the  floorplanning  stage. 

It  appears  that  80  %  -  90  %  of  the  computation  time  used  by  our  floorplanning  algorithm  is 
spent  solving  the  compaction  NLP.  Fig.  66  shows  the  average  change  in  objective  value  ( dotted 
line )  and  the  average  function  evaluation  time  ( dash  &  dotted  line )  needed  to  solve  the  compaction 
NLP  versus  the  number  of  subsystems  in  a  given  design.  We  created  our  test  cases  by  first  choosing 
a  common  set  of  typical  physical  property  values,  operating  conditions  and  performance  specifica¬ 
tions  for  each  subsystem.  Each  data  point  represents  the  average  of  20  randomly  generated 
instances.  One  standard  deviation  about  the  mean  is  also  shown. 

As  expected,  the  standard  deviations  of  both  time  and  objective  value  increase  as  the  number 

subsystems  increase.  This  is  because  the  solution  space  grows  in  proportion  to  (|N|!)2  •  8^  . 

It  appears  from  the  tests  we  have  conducted  that  our  algorithm  scales  reasonably.  This  can 
be  attributed  to  the  success  of  our  initialization  procedure,  which  helps  the  NLP  solver  converge  to 
locally  optimal  solutions  in  reasonable  time.  While  the  compaction  NLP  appears  to  scale  accept¬ 
ably,  the  SA  must  evaluate  the  NLP  more  often  as  designs  become  larger-scale  to  achieve  satisfac¬ 
tory  results.  We  have  produced  multiplexed  serpentine  designs  containing  up  to  30  subsystems  in 
under  one  day.  As  far  as  we  know,  these  designs  are  significantly  larger  than  any  multiplexed  ser- 
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FIGURE  67.  Example  routing:  (a)  Minimum  length  only,  (b)  Minimum  length  and  number  of  bends. 

pentine  designs  presented  in  the  literature.  However,  it  is  apparent  that  as  problem  size  continues  to 
increase,  the  computation  time  required  by  our  algorithm  will  become  prohibitive.  We  are  currently 
investigating  methods  that  will  allow  us  to  solve  larger  scale  problems  within  a  reasonable  amount 
of  time. 

MicroChannel  Routing  Results.  Fig.  67(a)  shows  a  minimum  length  routing  for  a  simple  5  sub¬ 
system  floorplan.  The  thick  lines  are  the  auxiliary  channels  connecting  ports  to  wells.  This  floorplan 
has  been  expanded  two  times  (i.e.  y  =  2).  This  means  that  there  are  at  least  three  possible  routing 
paths  between  each  subsystem  (each  subsystem  edge  and  one  path  in  between). 

Notice  that  since  there  are  many  degenerate  equal-length  paths  between  ports  and  wells,  the 
routing  solution  contains  many  bends.  Fig.  67(b)  shows  a  routing  solution  for  the  same  example 
when  both  the  channel  length  and  the  number  of  bends  have  been  minimized.  Typically,  both  the 
total  channel  length  and  the  port-well  assignment  remain  unaltered  between  the  minimum  length 
routing  and  the  minimum  length  and  bend  routing  (although  this  behavior  is  not  guaranteed).  This  is 
because  bend  minimization  represents  only  a  small  fraction  of  the  overall  objective  value.  In  gen¬ 
eral,  it  does  not  cause  re-routing  to  occur. 

While  bend  minimization  helps  to  reduce  the  number  of  non-unique  solutions,  it  increases 
the  optimality  gap  and  makes  the  problem  more  computationally  difficult.  In  our  experiments,  a 
standard  IP  solver  [117]  can  typically  prove  a  proposed  routing  grid  infeasible  in  less  than  30  sec¬ 
onds.  Therefore  the  decision  to  expand  the  grid  is  made  efficiently.  When  bend  minimization  is 
added  to  our  problem,  obtaining  a  proven  optimal  solution  often  becomes  computationally  prohibi¬ 
tive.  Furthennore,  no  routing  computation-time  correlation  can  be  drawn  based  on  the  number  of 
subsystems  in  an  arrangement  since  the  routing  problem  scales  with  the  number  of  grid  edges  and 
not  the  number  of  subsystems.  Currently,  we  cope  with  these  issues  by  allowing  routing  solutions  to 
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FIGURE  68.  Total  routing  length  vs.  number  expansions. 


have  a  small  relative  optimality  gap.  Despite  these  complications,  we  have  obtained  reasonable 
quality  solutions  in  1  hour  of  CPU  time  [118]  for  10  subsystem  designs. 

Fig.  68  shows  the  typical  operation  of  our  routing  algorithm  for  the  5  subsystem  example.  It 
illustrates  the  trade-off  between  minimum  channel  length  and  minimum  design  area.  Fig.  68  depicts 
how  iteratively  expanding  a  design  influences  the  total  auxiliary  channel  length.  For  expansion 
numbers  of  y  =  0  and  y  =  1,  no  feasible  routing  can  be  found.  A  feasible  routing  is  obtained  at  y  =  2. 

Normally,  we  would  terminate  the  algorithm  after  optimizing  the  first  feasible  routing. 
However,  if  the  design  is  expanded  a  third,  fourth  and  fifth  time,  we  notice  that  we  achieve  what 
appears  to  be  a  globally  minimal  length  routing  at  y  =  3.  This  occurs  because  for  low  y's,  the  routing 
grid  is  highly  constrained.  Therefore,  the  first  feasible  routing  solution  will  often  contain  long  chan¬ 
nels.  As  y  increases,  more  possible  routing  paths  are  created  and  shorter  routes  become  available. 
Eventually,  continuing  to  increase  y  will  result  in  longer  channels  simply  because  the  design  is  now 
larger  and  the  channels  must  traverse  a  longer  distance  to  connect  ports  to  wells. 

The  minimum  length  solution  at  y  =  3  is  also  significant  because  this  solution  most  closely 
approaches  the  solution  obtained  by  the  well  placement  algorithm  described  earlier.  Recall  that  in 
the  well  placement  algorithm,  route  continuity  is  enforced,  but  route  planarity  is  not  enforced.  The 
well  placement  algorithm  is  in  fact  a  lower  bound  relaxation  of  our  IP  routing  formulation.  We  are 
currently  investigating  how  to  utilize  this  infonnation  to  help  us  more  quickly  discover  optimal  or 
near  optimal  routing  solutions. 
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4.4.6.  Summary 

We  have  presented  a  design  automation  methodology  that  is  capable  of  generating  full-cus¬ 
tom  multiplexed  LoC's  for  CE  applications  in  only  minutes  to  hours.  We  have  adapted  and  com¬ 
bined  VLSI  circuit  design  techniques  with  optimal  design  methods  to  perfonn  LoC  design  and 
layout.  We  are  currently  investigating  methods  for  parallelizing  our  floorplanning  algorithm  to  han¬ 
dle  larger  problem  instances.  We  are  also  investigating  ways  to  reduce  the  number  of  edges  in  rout¬ 
ing  grids  while  still  maintaining  a  high  level  of  connectivity.  We  believe  that  our  experience  with 
multiplexed  design  provides  us  with  a  tool  set  that  can  be  extended  to  handle  multifunction  chips 
that  contain  many  complex  chemical  operations. 

5.  Conclusions 

The  primary  conclusion  of  this  work  is  that  it  is  possible  to  develop  biofluidic  microsystem 
models  for  optimization-based  design  synthesis.  The  separation,  mixing  and  reaction  models  devel¬ 
oped  in  this  project  are  capable  to  meet  the  needs  of  iterative  simulations  without  any  additional 
augmenting  infonnation  from  continuum  simulation.  At  the  start  of  this  project  the  state  of  the  art 
was  parameterized  lumped  element  models  for  simple  geometries,  and  reduced  order  models  for 
more  complicated  geometries.  This  project  has  extended  the  boundary  for  which  parameterized 
lumped  element  models  are  available  as  well  as  extended  the  applicable  range  of  design  dimen¬ 
sions.  For  example,  the  models  and  the  accompanying  simulation  capability  were  able  to  predict  the 
performance  of  a  highly  convective  dispersion  in  a  hybrid  electrophoresis  microchip. 

Component  and  system  optimization  using  these  models  were  demonstrated.  They  invari¬ 
ably  improved  the  quality  of  the  design  compared  to  the  existing  design  published  in  the  open  liter¬ 
ature. 

6.  Recommendations 

Further  development  of  application-specific  biofluidic  microsystem  design  should  proceed, 
the  suite  of  modeling,  simulation  and  optimization  tools  developed  in  this  project  are  now  a  basic 
design  framework  for  design  improvement.  An  NHGRI  supported  project  on  desinging  DNA 
sequencing  microchips  will  be  the  first  area  where  this  project’s  investigators  will  look  at.  Other 
DoD  targets  may  be  potentially  identified  by  one  of  the  graduated  students  from  the  project  who  is 
now  working  with  Northrop  Grumman  or  by  another  of  the  graduated  students  working  at  CFDRC. 
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